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ABSTRACT 

We study spherical and disk clusters in a near-Keplerian potential of galactic 
centers or massive black holes. In such a potential orbit precession is commonly ret- 
rograde, i.e. direction of the orbit precession is opposite to the orbital motion. It is 
assumed that stellar systems consist of nearly radial orbits. We show that if there is a 
loss cone at low angular momentum (e.g., due to consumption of stars by a black hole), 
an instability similar to loss-cone instability in plasma may occur. The gravitational 
loss-cone instability is expected to enhance black hole feeding rates. For spherical sys- 
tems, the instability is possible for the number of spherical harmonics ^ ^ 3. If there 
is some amount of counter-rotating stars in flattened systems, they generally exhibit 
the instability independently of azimuthal number m. The results are compared with 
those obtained recently by Tremaine for distribution functions monotonically increas- 
ing with angular momentum. 

The analysis is based on simple characteristic equations describing small pertur- 
bations in a disk or a sphere of stellar orbits highly elongated in radius. These char- 
acteristic equations are derived from the linearized Vlasov equations (combining the 
collisionless Boltzmann kinetic equation and the Poisson equation), using the action - 
angle variables. We use two techniques for analyzing the characteristic equations: the 
first one is based on preliminary finding of neutral modes, and the second one em- 
ploys a counterpart of the plasma Penrose - Nyquist criterion for disk and spherical 
gravitational systems. 

Key words: Galaxy: centre, galaxies: kinematics and dynamics. 



1 INTRODUCTION 

Mechanisms of "fuel" supply for galactic nuclear activity 
usually assume exposure of stars and gas clouds to a non- 
axisymmetric gravitational potential. The bar mode insta- 
bility and tidal action from nearby galaxies are most com- 
monly c onsidered to be respon sible for formation of the po- 
tential (|Sulentic fc Keellll990D . We believe, however, that 
the nuclear activity results from processes in the immedi- 
ate vicinity of central objects. It is unlikely that large-scale 
instabilities, such as the global bar mode, can provide for 
precise targeting of the star or gas flow towards the center. 

As an example of local mechanism that can maintain 
the activity we consider an instability in the stellar envi- 
ronment of the galactic center. This instability has a well- 
known prototype in plasma physics: the loss-cone instabil- 
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ity in the simplest plasma traps s imilar to mirror machines 
(see the pione e r pap er by Roscnbl uth. Pos3 |l965), and, e.g., 
iMikhailovskvl (|l974l )l. which is due to peculiar anisotropy in 
the velocity distribution of plasma particles. The anisotropy 
is caused by departure of particles with sufficiently small 
velocity component transverse to the symmetry axis of the 
system. The presence of this "loss cone" produces deforma- 
tion of the plasma distribution function (DF) in transverse 
velocities, giving it unstable (beam-like) character. 

Similar deformation in the DF in angular momentum 
takes place in clusters in case of deficiency of stars with low 
angular momentum due to their absorption by the galactic 
nucleus, black hole, or to other reasons. Then the deforma- 
tion can have a "beam-like" character: the DF becomes an 
increasing function of angular momentum, dfo/dL > 0. In 
this case, provided an additional condition discussed later 
is met, the deformation can trigger the instability which we 
shall call the gravitational loss-cone instability. 

We have mentioned the principal possibility of the 
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instability in, e.g.. iPol vachenk o fc ShukhmarJ (|l980l ). and 
iFridman fc Polyachenk o (1984). However, the search of a 
specific example of the gravitational loss-co ne instability 
has be en unsuccessful until one of the authors (|Polvachenkd 
Il991bl ) found the desired instability in a disk model. The 
delay and initial difficulty in finding it might be due to the 
unusual character of this instability, which is directly related 
to relatively very slow precession motion of stellar orbits. 
Thus typical frequencies and growth rates here are anoma- 
lously small, if measured in, for example, orbital frequencies 
of stars. Yet it does not mean slowness in absolute units tak- 
ing into account swift growth of all typical frequencies when 
goin g from the galactic periphery to its center. 

iTremaind (|2005l ) have studied the secular instability, 
which is identical to the gravitational loss-cone instability. 
He examined the disk and sphere models of a low-mass stel- 
lar system surrounding a massive central object. In such 
"near-Keplerian" systems, the gravitational force is domi- 
nated by a central point mass. For the disk models with ar- 
bitrary orbits, Tremaine has found unstable solutions. Note 
that for disks of nearl y radia l orbits the instability was 
proved by iPolvachenkol (|l991bh . and we will give here an- 
other p roof based on the gen eral integral equation for eigen 
modes (|E. Polvachenkdl2005l ) . Stability of spherical s ystem s 
for ar bitrary orbits has also been probed into by Tremaing 
l|2005l ). who found no evidence of instability for I ^ 2 modes. 
The study for I ^ 3 modes in general case is difficult, but 
it becomes feasible if one restricts consideration to nearly 
radial orbits. In this paper we show that the loss-cone insta- 
bility occurs just from 1 — 3 (see Sec. 4 and Sec. 5). 

For massive black holes in galactic centers, the coUi- 
sional diffusion and subsequent partial absorption of near- 
est stars is most often consid ered as a mechanism pro - 
viding for nuclear a ctivit y (e.g., iLightman. Shapirol l)l977h : 
IShapiro. MarchantI (|l978l ll. The very existence of the coUi- 
sionless (collective) mechanism may initiate revision of the 
dominating viewpoint regarding the nature of the activity. 
In this paper, however, we present a mere demonstration of 
the existence of gravitational loss-cone instability in simplest 
models. An exception is some preliminary estimations of ef- 
ficiency of the proposed collective mechanism in the Sec. 5. 

Existence of the above-mentioned additional condition 
for the instability originates from fundamental distinctions 
between gravitating and plasma systems. In gravitating sys- 
tems, particles have only one kind of "charge", and they 
attract each other. This fact ultimately leads to the Jeans 
instability substituting Langmuir oscillations in plasma (e.g. 
[Fridman &i Polyachenko 1984) . In systems with nearly ra- 
dial orbits we are going to study, there is a specific form 
of the Jeans instability ca l led the radial orbit instability 
|Polyachenko fc Shukhmax] Il98ll : I Fridman fc Polvachenkd 
1 19841 ). It develops only in systems with prograde orbit pre- 
cession (see Fig.[l}. Conversely, the gravitational loss-cone 
instability can occur only when orbit precession is retro- 
grade. This retrograde precession is the additional condition 
of the instability. 

As is well-known, the radial orbit instability is sup- 
pressed if the disper sion of orbit preces sion velocity exceeds 
some critical value () Polvachenkd [l99^ 1 . We shall obtain a 
similar result for the gravitational loss-cone instability (see 
Sec. 4). 

Below we study the gravitational loss-cone instabil- 



ity in two models representing active stellar subsystems 
with nearly radial orbits. A s noted abov e , the in s tability in 
disks was studied earlier in IPolvachenkol (|l991bh : [lYemainel 
(2005. ) . Nevertheless, we believe that it merits more detailed 
consideration here, although the main goal of this paper 
is to st udy instability in spherical systems. In IPolvachenkd 
(|l99ld ). instability has been determine d using the P en- 
rose - Nyquist general criterion (see, e.g., IPenrosd l|l960l ) or 
iMikhailovsk"^ (|l974l 'l). 

In this paper we apply a more illustrative method based 
on finding neutral modes and subsequent application of the 
perturbation theory. This allows us to obtain the frequencies 
and growth rates of perturbations corresponding to small de- 
viations from the neutral modes into the unstable domain. 
For unstable modes remote from the neutral ones, their com- 
plex frequencies are found numerically. At first we shall carry 
out this instability analysis for a simpler disk model (Sec. 3), 
then we shall use it for more complicated spherical geometry 
(Sec. 4). 

Another reason for re visiting disks is that a principal 
characteristic equation in IPolvachenkol l|l99 1bil was taken 
ready-made, without derivation. Here we present a detailed 
derivation of this equation. In addition, we justify the use 
of a suitable rotating spoke approximation: a spoke consists 
of stars with fixed energy E — Eq and low values of angular 
momentum L. The approximation is then applied to the 
spherical model. 

The paper is organized as follows. The characteristic 
equations derived in Sec. 3.1 and 4.1-4.4 are then applied 
to studying the gravitational loss-cone instability of disks 
(Sections 3.2 and 3.3) and spheres (Sections 4.5-4.8). In 
particular, in Sec. 4.8 we obtain general instability criterion 
for spherical systems analogous to the well-known Penrose - 
Nyquist criterion for plasma, and establish its correspon- 
dence to our neutral-mode approach. The instability of var- 
ious DF functions is discussed in terms of this criterion. 
The study is prefaced with an overview of the orbit pre- 
cession in the axial and centrally-symmetric gravitational 
fields (Sec. 2). Appendix A is devoted to derivation of a ba- 
sic integral equation for spherical systems in terms of the 
action-angle formalism, and in Appendix B we prove the 
instability criterion theorem. 



2 SOME REMARKS ON THE ORBIT 
PRECESSION 

In low-frequency perturbations of stellar clusters we are 
interested in, with typical frequencies, ui, of order of the 
mean precession velocity of near-radial orbits, these lat- 
ter participate as a whole (in contrast to high-frequency 
perturbations, for which lj is of order of orbital fre- 
quencies; they depend on individual stars). A detailed 
justification of these statements (w hich are fair ly obv i 
ous) can be found i n our papers |P olvachenkd ('l99! 
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and in lLTOden-BelTll979l ). For perturbations of interest, 
precessing orbits replace individual stars. So a preliminary 
overview of some useful data on precession becomes very 
desirable. 

In spherical pote ntials, star orbits are rose ttes that gen- 
erally are not closed l|Landau fc Lifshitj[l975 ) (see Fig.di). 
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Figure 1. Diagram explaining the physical mechanism of radial orbit instability: a - two typical radial orbits Aq and -Bo in equilibrium 
state, A and B are their positions in perturbed state from both sides of the perturbed potential minimum (bold dashes) at initial instant 
of evolution; b - stars at the orbits A (now A') and B (now B') have gained small angular momentum (L < at A' and L > at S'); 
c — a case of prograde orbit precession (the direction of precession is indicated by arrows) : further merging of large axes of orbits takes 
place, i. e. radial orbit instability develops; d - retrograde precession: return of large axes to the initial equilibrium position, i. e. neutral 
oscillations. These oscillations can be unstable in the presence of loss cone. 




Figure 2. Typical stellar orbit in the plane of axially-symmetric 
disk or in the spherical cluster a - "rosette trajectory" in the 
inertial system of reference b - closed precessing orbit in a system 
of reference rotating with angular velocity Q = Qpr. 



It is possible, however, to find a rotating (with angular ve- 
locity, n) reference frame in which the orbit is a closed oval. 
Therefore, star movement along the rosette are quick os- 
cillations in a closed oval, which in turn slowly rotates (or 
precesses) with the rate fipr — Q (Fig.[2j3). The latter is the 
orbit precession rate. 

There are only two potentials (| Arnold! Il989l ) in which 
any orbits are closed ellipses: (i) for quadratic potential, 
^o{r) oc r^, the ellipses are symmetric with respect to 
the center, and (ii) for the point mass potential, $o(?') = 
—GMc/r, in which the orbits (Keplerian ellipses) are asym- 
metric. The ellipses in these two potentials are resonance or- 
bits, since the ratio of the radial oscillation frequency, to 
that of the azimuthal oscillations, Q2, is 2:1 for the quadratic 
potential, and 1:1 for the Keplerian one. As the orbits are 
closed, the precession rate is zero. 

Small deviation from these particular potentials results 
in a slow precession with the frequency flpr much smaller 
than the typical frequencies of orbital motion, fli and . It 
occurs, for example, in the centers of galaxies. In the absence 
of a central point mass, the potential is almost quadratic, 
so the ratio of radial and azimuthal frequencies is close to 



2:1. If a point mass is present, this ratio is close to 1:1 in 
the central region. The deviation from the exact resonance 
(and thus slow precession) is caused here by gravity of stars 
around the central point mass. 

In smooth potentials, nearly radial orbits, with angu- 
lar momenta small compared, for example, to the angular 
momenta of circular orbits with the same energies E), are 
almost resonant for any amplitude of radial oscillations. This 
statement can be easily proved. Indeed, the azimuthal fre- 
quency, Q2, to radial frequency, ili, ratio for a star in the 
gravitational potential <l>(r) is, by definition: 

n2{E,L) _ Aip 



ni{E,L) n 



where 



dr 



Vr(r; E, L) 



dr 



^2E -2^r) - L'^/r 
(2.1) 



is the rotation angle of the stellar trajectory as its radius 
changes from rmax to rmin, E is the energy, and L is the 
modulus of angular momentum. 

Now let us calculate the asymptotic behavior of (|2.1|l 
for L ^ 0. To do this, we shall analyze the expression for 
radial velocity, — 2E — 2"I>(r) — /r^; its zeros define the 
turning points. 

In the case of a non-singular potential, $(0) is finite, 
and one may assume $(0) = 0. Obviously, the left turning 
point is rmin ~ L/\/2E. Since at low L the value of rmin is 
small, the main contribution to the integral comes from a 
region near the lower limit. So we have 



/■■ 



dr 



^2E - 2<i>{r) - /r^ 



dr 



2Sr2 
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Substituting {^J2E/L) r — x, we obtain 
L dr 



f 



dx 



and finally A(/p = 7r/2. Thus, 2/S.ip — n, and the orbit is a 
straight line which passes (almost) through the center. The 
absolute value of the ratio ^2/^1 is 1:2. 

Since for the Keplerian potential this ratio is 1:1, the 
question arises: Is there a continuous transition from the 
non-singular case, where the ratio is 1:2, to the Keplerian 
case? If such a transition exists, i.e. if the angle Aip can vary 
smoothly from 7r/2 to vr, a full stellar trajectory cannot be 
a straight line: a radial direction from the apogee to the 
perigee should change as a star goes from the perigee to the 
apogeeQ The angle of rotation is tt < 2 A(p < 2n, and the 
trajectory looks like spokes of a bicycle wheel. The number 
of spokes is finite if Aip/n is rational, otherwise the spokes 
will fill up a circle. 

Now let us consider, for the neighborhood of the center, 

$(r) = -^, 0<s<2. (2.2) 



The family of potentials (|2.2|l includes both the non-singular 
potential, corresponding to s = 0, and the Kelperian poten- 
tial corresponding to s = 1. Near the center, the absolute 
value of the radial velocity is 



^2S - 2$(r) - L ~ \/2ar-» - /r'^ , 

and thus rmin = (-£//2q)^'''^~°\ As before, the main contri- 
bution to the integral comes from the lower limit, so 

L dr L dr 



y^2£-2$(r) - LVr^ 



L2 



Substitution (2Q/L^)r^~° = x^ gives A<p = 7r/(2 - s), 
which leads to the frequency ratio 



HaM = 1/(2 -s). 



(2.3) 



The relation (|2.3p connects the non-singular and Keplerian 
potentials. We would like to stress again that in this case 
the stellar trajectory has a sharp turn almost in the center. 

Fig-E] shows a schematic trajectory in the potential 
<I>(r) ~ —r~^^^ (s = 1/2). One can see the sharp turn of the 
orbit with the rotation angle 2Aip = 2n/{2 — s) = 47r/3 = 
240°. Numbers 1, . . . , 6 trace a star in the trajectory. The 
star moves counter-clockwise, in accordance with the posi- 
tive sign of the angular momentum, L > 0. The trajectory 
is closed since 1/(2 — s) is a rational number (2/3). 

Further we focus our attention on spherical systems 
with near-radial orbits in two special gravitational poten- 
tials: 1) for a singular near-Keplerian potential, and 2) for 
an arbitrary non-singular potential. We shall refer to 1:1- 
orbits in the former case, 2:l-orbits in the latter case. Recall 
that the most obvious difference between these orbits is re- 
vealed in the degenerate case of radial motion: a l:l-orbit 
turns into a ray travelling from the center, while a 2:l-orbit 
turns into a line segment, symmetric to the center. 




Figure 3. Trajectory of a star with a small angular positive mo- 
mentum L. Parameter s = 0.5 



It turns out (see Introduction) that the gravitational 
loss-cone instability occurs if the orbit precession is retro- 
grade. Thus it is useful to have expressions for the precession 
rate, fipr, for both types of the orbits at hand. 

l:l-orbits. In the case of a low-mass spherical cluster 
around the central mass Mc, one can write 



dr 



2E + 2 GMc/r - 2$G(r) - /r2 



and 

ill TT 



dr 



Aip 



r2 yj2E + 2GMc/r ~ 2$G(r) - L'^ /r 



where Atp is the rotation angle of a star (cf. ((2T|) in the tra- 
jectory between rmin and rmax. If there is no precession, i.e. 
the contribution of $G(r) to the total potential is negligibly 
small, the angle would be tt. It is these small deviations of 
the angle that lead to slow rotation of the elliptical orbit, or 
precession: 



D,pr = D.2 - 0,1 = Qi {A(fi/n - 1) 



(2.4) 



The expression for the precession velocity has obvious sense. 
After rewriting it in the form 



Op 



T/2 ' fii' 



one can see, that if during one half oscillation, from rmin to 
''max (for the time T/2), a star travels an angle exceeding tt, 
it would imply that the apogee drifts with angular velocity 

Our goal now is to derive an expression for precession veloc- 
ity of order 0{$g)- Since Aip/Ti— 1 = 0{^g), it is sufficient 
to retain the 0(1) order only in calculating the Sli factor in 
(P^ . Then we have 



A^^LA[(2li5l)i/2 / ^^^b-r){r-a) 



dr 



^ Note that in the Keplerian limit radial orbits degenerate into 
a "ray". More precisely, the incoming and outgoing rays merge 
together. 



(2li51)- 



-1/2 



<i>G('") dr 



r \J {b — r){r — a) 



(2.5) 
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Here we introduced new variables a{E,L) = rmin{E,L) 
and b{E,L) = ruia,x{E,L), instead of E and L, where 
a,b = GM,/{2\E\) T {[GM,/{2\E\)]' ~ Ly{2\E\)y^\ 
For later purposes, we need some useful relations valid in 
the Keplerian potential: 



, GM , 



L = 



2GMc at 
a + b 



(2.6) 



One can show that the first item in r.h.s. of (12. 5p is zero, so 
one obtains for the precession velocity 



n{E) d 
TV dE 



[l{2\E\)-''^ [ — 



<l?G(r) dr 



{b — r)(r — a) 



where Q,{E) = Q.i{E) = Q.2{E) ^ {2\E\f''^ / GM^. Using 
(|2.6|l . it is easy to obtain an expression for precession velocity 
in variables (a, b): 



a \ 

o 

[Ubf 



.2d_ 
da 



^G{r) dr 



i r^{b~r){r-a) 



(2.7) 



where Q. = (GM^)^''^ [2/(a + 6)]^''^ The expression for D,pr 
can be simplified, assuming a to be small (nevertheless, we 
retain a in the denominator of the argument of the square 
root): 



a 



»(b) ,2 d 
^ TT GMc db 



[vrbj 



^air) dr 



i r ^ [b ~ r)[r - a) 



where o}^"^ = L/ J 2GMc . It is possible to show that integral 



Jia,b) = I 



^G{r) dr 



i r^{b-r){r-a) 

b 

J{a,b) = I $G(r)d[-l= 



may be written as 



(b + a)r - 2abl 



r{b~a) J • 

a 

Integrating by parts and then differentiating over b yields: 



d_ 

db 



'ab J {a, b) = 



d^c 



dr b — a 



1^3/2 



Finally, we obtain for near-radial orbits 

^p,{E,L)^vj{E)L, 

2 



■ dr ■ 



d^G 
dr 



■ dr. 



vj{E) = - 



TxGMa b 



d^G 
dr 



b — r 



■ dr 



(2.8) 



with b = GMc/\E\. From Eq. (I2.8|l . it is clear that for such 
orbits and potentials the precession is always retrograde 
since ^^(r) — GMair) /r^ > 0, where MG{r) is the mass 
within a sphere of radius r. (Moreover, it is easy to show 
also from (|2.7|l that for near-Keplerian orbits precession is 



retr ograde eve n for the case of arbitrary eccentricity (see 
also lTrema"in3 l2005')). 

Thus, for spherical near-Keplerian systems, the retrograde 
precession is common. As we will see later, for 2:l-orbits, 
the situation is vice versa: in real potentials, the precession 
is prograde, whereas the retrograde precession occurs in sys- 
tems with exotic density distributions. 

2:l-orbits. There are several different expressions for the 
precession velocity of near-radial orbi ts of this type. The 
expression given m IPolvachenkol (Il992l ) is somewhat incon- 
venient for practical use, because it contains a procedure of 
passage to limit. We shall give a more convenient formula. 
For the precession velocity, one has: 

n.(i.,L)^.i(i.,L)[g^-i]. 

The goal is to calculate the derivative of the precession ve- 
locity at constant E and L — 0, i.e. [d^lpr{E, L)/dLj 
Taking into account that f2i(£',L) = Q{E) + 0{L'^), one 
obtains 



^lpAE,L)^n{E) (Av9/7r-i) 



where 



dr 



72£-2$o(7-)-LVr2 



is the rotation angle of a star in the trajectory between rmin 
and r-jnax, which is equal to i vr for any non-singular poten- 
tial, in the leading order in L. The next (linear) term in 
expansion of Aip gives a value of the precession velocity and 
defines the direction of precession (prograde or retrograde) , 
for a highly elongated orbit, with a small angular momen- 
tum. It is clear that if A(p > n/2 (or 2Alp > n), the preces- 
sion direction coincides with the rotation of a star (prograde 
precession), and vice versa. Thus, the sign of the derivative 
d/dL {Aip — n/2) determines the precession direction. Let 
us denote (j}(E, L) — Aip — 7r/2, so that 

zu{E) = [dQp,{E,L)/dL]^^^ = ^ [d^{E,L)/dL]^^^. 

(2.9) 



We have 



r,„ax{i5,L) 



Aip : 



dr 



r-minCB.-L) 



^2£-2$o(r) -LV^ 



r,„ax {i5,L) 

dE J r2 
Integrating by parts yields 



^2E -2<S>o{r) - L'^/r'^ dr. 



Aip = 



_d_ f 

dE J r* 



dr 



\/2£ -2$o (r)-LVr2 
d f "I>o(r) dr 



dE J r ^2£-2$o(r)-LVr2' 

One can show that the main contribution to the first in- 
tegral comes from the lower limit, thus one can extend the 
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integration to infinity. Since, for non-singular potentials, one 
can neglect ^o{r) ^ $(rmin) (as compared to E) in the root 
argument, we obtain 



I 



dr 



^2£-2$o(r)-_LVr2 



2E 



dx 



1 X' 



so for 4'{E, L) we get: 
d 



-L- 



dE 



dr 



r„ax(B) 

^ ^ n dE J 



r 




find 






dr 


r 


,/2£-2$o(r) 



(2.10) 



There is also an alternative expression, without differentia- 
tion with respect to E. Omitting details, we give the final 
expression: 



dr 

72 



1 



2E - 2$o(r 



'2E 



,iE)V2E 



(2.11) 



From the derived expressions it follows, for example, 
that for the expansion of a potential, $o(r) = $o(0) -I- 
Qq (^^ r'^ + P r'^ + ■ ■ ■) , at small r, the precession direction 
is determined by the sign of /3, since = — /3. Thus, stars 
with small radial amplitude have retrograde precession if (3 
is positive. From Poisson's equation one immediately has 



47rGr2 dr dr 



M(l + ^/3r^ 
47rG V 3 



i. e. density increases with radius. Such behavior is unreal- 
istic, and we conclude that the gravitational loss-cone insta- 
bility is impossible in spherical clusters with non-singular 
poten tials. Note that in d isk systems this instability is pos- 
sible (|Polvachenkdll991bD . 



3 GRAVITATIONAL LOSS-CONE 
INSTABILITY IN DISKS 

3.1 Derivation of an equation for eigenmode 
spectra 

The simplest relations for describing the instability under 
study can be obtained in a model with an active stellar sub- 
system in the form of a disk with nearly-radial orbits. The 
aim of this Section is to derive characteristic equations for 
small perturbations in such a model. The derivation involves 
a series of successive simplifications of the initial linearized 
Vlasov equations. 



Omitting a number of those steps, we take, as the start- 
ing point, the following integral equatior0: 

. G f f dE'dL' 



E 



2tt J J ni{E',L') 

[mfl2{E',L') +1' Q.i{E',L')] dF/dE' + mdE/dL' 

> 

''=-°° uj-mn2{E',L')-l'ni{E',L') 

xUij,iE,L-E'L')<i>t,iE',L'). (3.1) 

Here G is the gravitational constant, the frequencies Qi — 
dHo/dE, or more briefiy, 0(1) = dHo/dl, Ho (I) is the 
Hamiltonian of a star in unperturbed state, I — (Ji,/2) are 
the actions, lj is the frequency of a perturbation, E and 
L are the energy and the angular momentum, respectively, 
F = F{E, L) is the unperturbed DF, and the Fourier compo- 
nents <!>;(£', L) and the kernel Tliji{E,L;E',L') are defined 
as follows: 



^i{E,L) = — / dwi^r{E,L;wi)] cos[lwi+mxiE,L;wi)], 

where $(r) is the radial part of the perturbed potential 
■^{r,ip-t) = $(r)e"*"'+''"*', t is the time, ip is the polar 
angle, m is the azimuthal index, wi is the angle conjugate 
to the radial action Ii, I is the integer index. 



Ili,ii{E, L; E , L ) = / dwicos[lwi + mx{E,L;wi)] 

J —7T 

dw'i cos [ I'w'i + mx{E' , L'; w'l)] 



x^[r{E,L;wi),r'{E',L';w[)], 

, f^" ,„ cos me 

'tp{r,r)^ de 

Jo ,/r'^+r'^—2rr'cost 



The function x is defined by 

r 

x{E,L-wi) =n2{E,L) f - 

J Vr 



dx 



>-mi„{B,L) 



{E,L;x) 



dx 



x'^Vr{E, L; x) 

(where Vr is the radial velocity of a star), or x(E,L;wi) 
(fl2/fii) wi ~ Sip, where 



r 

5ip = L J 



dx 



x^Vr{E, L; x) ' 



'■mi„{B.i) 

For definiteness, we consider systems consisting exclu- 
sively of 2:l-orbits in this Section. However, the final form 
of the desired equation for the case of l:l-orbits (i.e., for 
near-Keplerian systems) is practically identical to that for 

2 Actually, 1 13.111 is a set of connected integral equations, but for 
brevity we shall call it simply 'equation' hereafte r. Il3.lt was ear- 
lier obtained in a paper by one of the authors (HI Polvachenkol 
using the actions — angles formalism suitable for problems 
of this type. 
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2:l-orbits. It is evident that slow modes, with angular rates 
of the order of the typical precession velocity, in the systems 
with 2:l-orbits are possible only when the azimuthal index 
m is even. In case of l:l-orbits the azimuthal number m is 
arbitrary. 

Although the equation (|3.1|) is an exact linear integral 
equation that allows one to determine the spectrum of eigen- 
modes for an arbitrary distribution of stars in the disk, for a 
stellar system with low angular momenta we are interested 
in, this equation is inconvenient as the system can include 
stars with both direct (L > 0) and inverse {L < 0) orbital 
rotation. The point is that the orbital frequency ^12{E, L) is 
discontinuous at L = 0: Q2{E,L = ±0) = ^ sl iii{E, L = 
0), sl = sign (L). So all items in the sum over I' are also 
discontinuous at L' = 0. The function itself and the 

kernel of the integral equation 11; ;/(L, L') are discontinuous 
too. (Hereafter we shall omit the arguments E and E' , pro- 
vided this creates no difficulties.) The discontinuity is very 
inconvenient. In fact it arises from a poor choice of the angle 
variables in the set of actions - angles, for the problems of 
interest involving stars with L < (note that discontinu- 
ity is quite i nessential in problems with nearly-circular or- 
bits (see, e.g.. lE. Polyachenkdl2005h . However, this difficulty 
is fictitious, and actually the equation can be transformed 
into a continuous form by means of a proper procedure that 
involves shifting of indices and transformation of the func- 
tions. The procedure is simple but cumbersome, so we give 
the final form of the integral equation without going into 
detail: 



n - ^ r f dE'dL' 
AE,L) - — 



E 



n' — — oo 



2-K j J Qi{E',L') 
{dF/dL')^g + (n'/m) ni{E' , L') dF/dE' 

Q.p - npr{E',L') - {n'/m) Q.i{E',L') 

xR„_^,{E,L;E',L')<t>n'{E',L'). (3.2) 

Here 0.^, dF/dE + dF/dL = {dF/dL)^g is the so called 
Lynden-Bell derivative of the DF which is by definition a 
derivative with respect to the an gular momentum at fixed 
Lynden-Bell's adiabatic invariant l| 19791 ) , Jj = Ii + ^ \I2\- 

(dF/dL)^^ = (dF/dL) = -i SL (dF/dh),^ + [dF/dh),^ 

and the precession velocity is 

9.^,{L')^Q.2{L')-\sL'ni{L'). 

Note that the function Jlpr (L) is continuous at L = (passes 
through zero). In the equation (|3.2p . the new function (con- 
tinuous at L = 0) 



K{E,L) = —j dw^^[r{E,L-wi)] 



and the new kernel continuous at L = and L' = 

dwi cos [null -I- mx{E, L\ wi)] 

-7r 

dw'i cos [n w'i~\- m x{E' , L'; w'l)] 

xi;[r{E,L;wi),r'{E',L';w[)] (3.4) 



dx 



appear. In (|3.3p and (|3.4[) 

r r 

x{E,L-w,)=n^,{E,L) f —p-—~L f — 

J Vr{E,L-x) J x^Vr(E, L;x) 

that is x{E, L;wi) = x{E,L; wi) - | sl wi. 

Then the equation (|3.2|l can be treated as the initial 
exact integral equation. Since below we shall concentrate on 
studying distributions localized in the vicinity of L = 0, it is 
important to keep in mind that the functions <j)n{E,L) and 
Rn, n' {E, L; E' , L'), and also the frequency fli{E, L) are not 
only continuous, but smooth at L = L' = as well: 

(f)n{E, L) = <!'„(£', Q) + a„{E)L + .. ., 
ni(£,i) = f2i(S,0) + O(L'), R„,n'{E,L-E',L') = 



, {E, 0; E', 0) + Pr.{E) L + f3„, {E') L' + , 



Although proof of this statement is rather non-trivial, 
we leave it beyond the scope of the article. Here it is par- 
ticularly important for us that the coefficients a„(i5) and 
Pn{E) for 11 = Q and n' = tend to zero, hence 

ME,L)=<t}Q{E,Q)+0{L''), 

Ro,o{E,L;E',L') = Ro,oiE,0; E' ,0) + 0{L') + OiL""). 

(3.5) 

It is precisely this fact that allows us to consider the kernel 
-Ro, and the function <j)o to be constant at the DF localiza- 
tion scales (with respect to L), when studying slow modes. 

Our following step is to transform (|3.2|l into an equation 
describing slow disk modes with frequencies of the order 
of precession velocities. The latter are always less than the 
orbital frequencies Qi and ^2, and for nearly-radial orbits 
f2pr <^ f2i, ^ 2. As it was exp l ained in considerable detail 
in a paper bv lE. Polvachenkg (|2004l ). for slow modes, only 
the items with n' = n = dominate in the sum of p.2|) . 
since they have minimal denominators. If these items are 
only taken into account, we obtain the equation 



ME,L) = ^ 



dE' dV 



2-K j j Qi{E',L' 
idF/dL')^^ 



np-Qpr{E',L') 



Ro.o{E,L;E',L')ME ,E). (3.6) 



X cos[nwi + mx{E,L;wi)], (3.3) 



This is the desired equation for slow modes. 

The next step consists in transforming (|3.6p into an 
equation convenient for a disk model with nearly-radial or- 
bits. Let us assume that in the domain of small angular 
momenta, the DF is F{E,L) = f '^\E),f^\L). The scale 
of localization domain, 5L, for the function f'^'"\L) near 
L = is assumed to be small. In general, the exact mean- 
ing of this smallness needs to be refined, but in any case 
the scale must be smaller than the characteristic length of 
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variation in momentum, AL, for all the functions appear- 
ing in the equation (|3.6p . i.e. (j)o{E,L), Ro,o{E, L; E' , L'), 
^li{E,L), flpr{E, L). The characteristic scale of variation 
for these functions is determined exclusively by the be- 
havior of unperturbed potential. As for the latter, we do 
not suggest any peculiar behavior and consider this po- 
tential to be non-singular at r = 0. So we can assume 
that 6L <C AL. In this case, due to the relations p.5|) . 
in the localization domain of the function /'■^■'(L), we can 
take ME,L) ^ <E>o(£,0) ee $(£;), Ro,oiE, L; E' , L') ^ 
Ro,o{E,0;E',0) = P{E,E'), Qi{E,L) ^ ni{E,0) = Q{E). 
For the precession velocity in a domain of small values of L, 
we have f2pr(S, L) « vj{E) L, vj = [dQ.p^{E, L)/dL] 

As a result, the two-dimensional integral equation (|3.6p 
reduces to a one- dimensional equation with the kernel de- 
pending on E and E' only: 



„ df ^\L')/dL' 



X / dU 



Hp - uj{E') L' ■ 



(3.7) 



To find the spectrum of eigenmodes, a numerical solution is 
required. However, one can predict immediately some qual- 
itative consequences. For P{E,E'), we obtain 

PiE,E')^mE)niE') /-^ l^^^^iry) 



with Vr{E, r) = sj2E-2^o(r) . It can easily be shown from 
(|3.7|l that in the "cold" case of purely-radial orbits, i.e. for 
/(^)(L)=5(L), 

It is easy to verify (see, e.g., IPolvachenkd [l993) that 
P{E,E') is a positive quantity. Then it is evident that in- 
stability or stability depends exclusively on the sign of i:i7(i?). 
Namely, instability occurs when it is positive, vj > 0. This 
is the radial orbit instability. Recall that near the center, 
where the potential <l>(r) = $(0) + r2o(| + 13 r'^ + . . .), the 
quantity zu is positive when /3 < 0, so that instability must 
occur. Note that such behavior of the potential is typical for 
most surface density distributions decreasing with radius. 

To investigate the spectrum of eigen oscillations in more 
detail, let us add another simplifying assumption. Namely, 
let us consider a model with monoenergetic distribution over 
energy, i.e. F{E, L) — S{E — Eo) f{L). In this case the inte- 
gral equation (|3.7|l reduces to a simple characteristic equa- 
tion for the complex (generally speaking) velocity fip: 



G 



2-kQ.{Eo 



■P{Eo,Eo) 



df{L)/dL 
— zu{Eo) L 



dL. 



One can turn from distribution over the angular momentum 
L to that over the precession velocities: fipr — ■oo{Eo)L = u. 
After denoting /(L) = / iv jvo) = fo{v), we obtain 



G 



2-kQ.{Eo) 



P{Eq,Eo) sign (zu) 



dfo{v)/du 

fin — V 



dv. 



Hereafter we are primarily interested in the case of retro- 
grade precession, w < 0, therefore let us write the charac- 
teristic equation in the following final form: 



-A 



f dfoju 
J ^P- 



)ldv 



dv, A — 



G 



2ttQ.{Eo 



or equivalently. 



l=A 



I 



{Vtp 



■ du. 



■P{Eo,Eo) > 0, 
(3.8) 



(3.9) 



It is easy to check that the equat ion (|3.9[l coincide jf| w ith the 
characteristic equation obtained IPolvachenkd ()l991b[ l in the 
so-called "spoke"-approximation. The derivation above is in 
effect a formal justi fication for this equation obtained earlier 
by V. Polyachenko (|l991bl ) using a semi-intuitive approach. 
In this approach, a set of stars moving along the same elon- 
gated orbit is regarded as a new elementary object replac- 
ing individual stars, and the dynamics of stars reduces to the 
dynamics of spokes (for slow processes) ; for an exten ded dis- 
cussion see the relevant papers by V. Polyachenko (|l991al : 
Il99ld ). The advantage of the spoke approach is that it is 
much simpler than the general methods commonly used. It is 
this approach that is appropriate for studying low-frequency 
oscillations and instabilities. However, its rigorous justifi- 
cation required the above, rather cumbersome calculations. 
This procedure was however necessary in order to make sure 
that the spoke approximation is reliable. These calculations 
were also useful in that they illustrated suggestions and as- 
sumptions required for the approach. 

Now we can proceed to the study of the resulting char- 
acteristic equation. 



3.2 Loss-cone instability of a disk in the spoke 
approximation 



Let us represent Eq. (3.8) in the form 



Q= du 



df(u)/du 



A'^ > 0, (3.10) 



where for a distribution with a loss cone, i.e. with a defi- 
ciency of stars with low angular momenta, the DF f{u) has 
a zero minimum at v = Q: 



/(0)=0, /'(0)=0, /"(0)>0. 



(3.11) 



We assume that this minimum is unique, and the DF looks 
like what is shown in Figs. 4 or 5. An important point is 
that the quantity Q in Eq. (|3.10|l is positive only when the 
orbit precession is retrograde. It is this circumstance that 
allows us to lean upon the analogy with plasma and use the 
formalism of the plasma theory when studying the insta- 
bility. Recall that for the case of direct precession, when 
Q < 0, Eq. (|3.1U|I describes only the radial orbit insta- 
bility (the modification of the Jeans instability for very 
elongated orbits) leading to the spokes mergin g together 
(|Polvachenko fc Shukhmadll972l : lAntonovlll973l). Based on 
the Penrose - Nyquist criterion (see, e.g.. |Penros3ll96d . or 



^ Exce pting an inessential slip in the coefficient A in lPolvachenkd 
lll99ld) . 
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Figure 4. Loss-cone instability in disk, a — symmetric distribu- 
tion over precession velocity 113.131 1 with a = 0; 6 - dependence of 
Re(Qp) and Im(f2p) on parameter Q. 



lMikhailovskvll974l ') . [Polvachenkol l|l991bl ) showed that a new 
instability must occur in systems of orbits with retrograde 
precession. Here we reproduce this proof using a somewhat 
different language which is more physically evident and even 
more constructive: e.g., it allows us to determine, in a rel- 
atively simple manner, the instability boundary in the pa- 
rameter Q. This approach is based on considering neutral 
modes. 

The essence of the approach is as follows. Suppose that, 
when the parameter Q changes, the initially stable system 
becomes unstable. This means that some value of Q = Qc 
exists for which a neutral mode appears. For such a neutral 
mode, the location of the resonance should coincide with an 
extremum of the DF f{v) - otherwise Eq. (|3.10p will have 
a pole on the real axis, bypassing which would necessarily 
result in an imaginary contribution into the right side of 
IpTTO)) : so that Eq. ((3!T0|l could not be fulfilled. For distri- 
bution with one minimum and two maxima, one can show 
that (i) a neutral mode with a resonance at the location of 
a higher maximum cannot exist; (ii) when one maximum is 
much higher than the other, and two maxima are sufficiently 
separated from each other, the neutral mode corresponding 
to the lower maximum is possible, this mode belonging to 



/ (x) =x^exp [-(x-0.5)^] 
0.8 

0.6 

0.4 





-3 



Im (i^) 0.4 




Figure 5. The same as in Fig. 4 for asymmetric distribution 
l|3rT3)l with a = 0.5. 



the same oscillation branch as the neutral mode related to 
the minimum of the DF (see Sec. 3.3). Consequently, for the 
neutral mode at the minimum, Q.p = G, while for the other 
neutral mode, Q.p = u\, where = u\ \8 the location of the 
lower maximum. 

Putting f2p = in (|3.10[) . we find the value of Q = 
for which the neutral mode connected with the min- 
imum can exist: 



Qc = 



(3.12) 



Here we have no need to elucidate the meaning of these 
integrals. Due to the condition (|3.1ip . they converge at v — 
in the ordinary sense. Obviously, the right side of (|3.12|l is 
positive. So Qc™'"' always exists, for arbitrary locations and 
heights of maxima. As to the sign of Qc™''^'' , it can be either 
(as mentioned above). If it is positive, the second neutral 
mode exists. Note in addition that one more neutral mode, 
with Qc ~ always exists; formally, it corresponds to the 
resonance at infinity. 

With a knowledge of the values of Qc for neutral modes, 
one can determine domains of instability in the parameter 
Q {Q > 0) as Qc are the margin values. For this purpose, 
let us apply the perturbation theory. 

First we consider the region near the boundary Q — 



Qc"^ . Let us defied from critical value: Q 



■){min) 



+ SQ. 
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The frequency Sip acquires an addition S^p, which is sim- 
ply equal to Q.p. We would like to find out what direc- 
tion should be taken to reach instability. Varying the quan- 
tity Q relative to Qc"""' and Q,p relative to zero in (l3.1Up : 
5Q = SQp du f' (u). Here the pole in the integral 
must be bypassed below. This gives 



SQ = 



oo 

SQp [jdvv-'^f{v) + iTTf'{Q)\. 



OO 

Denoting A = -j- dvu~^ fi^)^ ^ = '^/''(O) > O5 obtain 
—00 

Qp = 5Q.p = SQ/{A + iB) = [{A-iB)/{A^ + B^)] SQ. 

As -B > 0, the instability appears when Q decreases below 
the critical value 

For the second neutral mode (with the resonance at 
the location of lower maximum), we use the same proce- 
dure of the perturbation theory to find that the instability 
is possible when Q is above the critical value Qc™'"'', as 
B = 7r/"(«i) < in the maximum. 

As a result, we conclude that in the absence of the 
neutral mode associated with the lower maximum, i.e. if 
Q(max) ^ unstable domain lies in the range < Q < 

Q(min) gygjj a, neutral mode exists, the range of instabil- 
ity becomes Q*™'"'' < Q < Qi™'"'. (It can be shown that 
Q^^^) < Q 

As for the question of where the resonance shifts when 
the parameter Q is deflected from the corresponding margin 
value into the unstable domain, it depends on the sign of A, 
defined by the integrals in the sense of principal value. In 
principle, these integrals can have any sign. In the particular 
case of symmetric DF, showed in Fig. [J] the quantity A is 



zero (for Q deviating into the domain Q < Q 



(min)N 



Stability 



(or instability) is governed only by the sign of B, and this 
does not depend on the sign of A, which can be associated 
with the angular momentum of the wave. 

We would like to say in this connection that here it is im- 
possible to consider our instability in terms of exchanged an- 
gular momentum between the wave and the resonance stars, 
as is done in plasma physics or in the theory of galactic struc- 
tures which has undergone appreciable developmen t sinc e 
the well-known paper by Lynden-Bell and Kalnajs (| 19721 ). 
The language for explaining the instability uses considera- 
tions operating with the momentum exchange between the 
wave and resonance stars. For instance, if the wave momen- 
tum is positive and the resonance stars lose their momentum 
transferring it to the wave, the momentum of the latter in- 
creases. This is the instability. However, such a language 
does not work in the case under consideration. The reason 
is that our systems are not weakly-dissipative, as is usual in 
plasma. In the latter we have a well-defined wave. All stars 
contribute into the wave dispersion properties, while the dis- 
sipation is determined by a small portion of resonant stars. 
Consequently, the dissipation is only a small correction. But 
now we have a completely different situation: the dissipation 
and dispersion parts of the wave (i.e., roughly speaking, the 
imaginary and real parts of dielectric permittivity) are of 
the same order. So our instability is not kinetic, in the ordi- 
nary sense, and such considerations do not work. However, 
in the case when the maxima are sufficiently separated, and 



their heights strongly difi^er from each other, we return to the 
usual, weakly-dissipative situation. Then the sign of phase 
velocity (i.e., the sign of A) must correlate with the incli- 
nation of the DF. In the following Sec. 3.3, we consider this 
case as well. 



3.3 Neutral modes with resonance at DF maxima. 
Investigation of stability in a model 
two-humped distribution 

To illustrate the above reasoning, we shall study a model 
example with DF 

fa(iy) = {u^ /vt) exp [-{u - aur)^ /ut] 

= 2:^ exp[— (x — a)^], (3.13) 

where we can control the maximum locations and heights. 

But first we shall prove the statement formulated in the 
preceding subsection that Qc for the neutral mode related 
to the higher maximum is always negative, while Qc for 
the lower maximum can have any sign. Suppose that the 
distribution has one maximum, similar to that showed in 
Fig. 4. Let the left maximum he a.t v = u\, and the right 
maximum at v = U2, the right maximum being larger than 
the left one: f(u\) < /(M2). We shall consider these two 
variants separately. 

(i) First let us assume that the neutral mode has fre- 
quency Sip — U2. Let us rewrite (|3.1Up in the form 



U2 



dp 



[f{^)-f{u2)] ^ 
{u~U2Y 



(3.14) 



Note that all integrals here can be considered in the usual 
sense since no problems arise concerning their convergence 
at V = U2. Since for the higher maximum, v = U2, /(f) < 
f(u2) everywhere, then the right side of (|3.14|l is obviously 
negative. Thus, we proved that there cannot be a neutral 
mode related to the higher maximum. 

(ii) Let us now suppose that Sip = u\. Let us split the 
integral into two parts: 



l(max) _ / rf/('^)/rft^ ^ Z"^^, 



V — U\ 



V — U\ 



{v - UiY 



Obviously, the integrals in the right side have opposite signs, 
so that the resulting sign can be either. In the case of the 
model p.l3p . the quantity Q^™^''' can be readily calculated: 



Q^r^'Ha) = 2^ [a^ia^ + l-i {a' + 1)], 

so that at a > flc = 2~^^'^ ~ 0.71, the quantity Q^^'^^\a) 
becomes positive. Note that in the given model, Q^"""' does 
not depend on a and equals y^yr. 

Fig-Slshows the dimensionless complex phase velocity of 
the wave. Sip, as a function of Q (in units of vt, for the model 
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Figure 6. Unstable region on Q-a - plane for the model l|3.13|l is 
bounded from above by the straight line Q = Q^™'"^ = from 
below by the curve Qi"'"''(a) = [a^\a'2 + 1 - \ {o? + 1)] 

(shaded) . 
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Figure 7. The tracks of eigenvalue xq in complex XQ-plane for 
various values a, indicated near corresponding curves. 



(|3.13[) with a = 0. It is evident that owing to symmetry the 
real part of the frequency Q.p (and the quantity A as well) 
is equal to zero. 

Fig.[5] shows Q,p{Q) for the model ((3TT3|) with a = 0.5. 
Since a < a^, then this model (as is the case with the model 
with a — 0) sees only a neutral mode related to the mini- 
mum. As the model is asymmetric, the real part of the fre- 
quency is not zero. It is negative for \SQ\ <C 1 as A > 0. 
Calculation shows that it remains negative when Q is far 
from the instability boundary. 

For sufficiently large values of a (a > Oc = l/\/2), the 
DF maxima (|3.13|l will be highly separated, and one of them 
becomes much higher than the other. Then qI™'"'', for the 
neutral mode corresponding to the lower maximum, becomes 
positive, and a second neutral mode appears. As a result, the 
instability domain looks like a horizontal band (converging 
with increasing a), between these two neutral modes (see 
Fig. 6). Fig. 7 shows the tracks of a complex eigenvalue xq 
in a complex xo-plane for various values of a. When the 



parameter Q decreases from Q 



(mill) 



r^/^ to Q('"""'(a) the 



position of the point on the corresponding curve changes so 
that Re(2;o) moves from to 2;inax(ffl), where a;niax(a) is the 
position of the lower maximum, while Im(xo) > and tends 
to zero on both ends of the curve. 

Based on the plasma analogy, there is no difficulty in 
understanding the physical essence of the instability under 
sufficiently large values of a. When a ^ 1, the DF becomes 
identical to the DF of plasma particles with a weak beam 
moving at a rate significantly higher than the thermal ve- 
locity of particles in the main plasma (a beam at the tail). 
Then the instability degenerates into the well-known beam 
instability. It occurs when the wave phase velocity is on the 
slope of the beam DF oriented towards the main plasma. For 
our model (|3.13|l with a > the phase velocities of unstable 
modes must be negative, which is the case in our calculations 
(see Fig. 7). 



4 LOSS-CONE INSTABILITY IN 

SPHERICALLY-SYMMETRIC SYSTEMS 

4.1 Basic set of integral equations 

As we did in the case of disks, here we start with an ex- 
act equation governing perturbations in a spherical system 
with DF F = F{E, L), where E = | + f L^/r^ + 'i'oir) 
is the star energy, L — rvj_ = r {vj + v^)^^^ is the ab- 
solute value of angular momentum, 'I>o('') is the unper- 
turbed gravitational potential. We assume that the DF 
does not depend on the third integral Lz = rii^sin6'. 
As is well-known, the spectrum of eigenvalues ui in this 
case is independent of the azimuthal number m. Thus, in- 
stead of a general representation of the potential and den- 
sity in the form of the sectorial harmonic $(t; r, 9, ip) = 
x{r)Yr{eA)e-'^' and p{t-r,e,^) = p{r)Yr{e e'^-\ 
we can restrict our consideration to a simpler vari- 
ant ^{t;r,e,p) = x('^)-P;(cos6')e~"*", p{t-r,e,ip) = 
p{r) Pi{cosd) e~*'^', where Pi{x) is the Legendre polynomial. 
The derivation of the basic integral equation (a set of inte- 
gral equations, to be precise) is also based on the action- 
angle formalism (as in the disk case) and presented in Ap- 
pendix A. Note that this formalism was fi rst used in the 
paper by Polyachenko & Shukhman (|l98ll ) as applied to 
spherical gravitating systems. 

Thus, the initial exact set of integral equations has the 
form {hjl'i = — oo, . . . , oo; ^2,^2 = —I, ■ ■ ■ ,1)'- 



Xii,i2iE,L) = 



4nG ^ ^ li, f dE' LdL' 
21 + 1 ' J Qi{E',L') 



xxi[i'^{E',L')nt^j^.^i,^^i,^{E,L;E',L')x 
^l[ni{E',L') +1'2 n2(-E",L')] dF/dE' dF/dV 



-l[ni{E',L')-li, Q.2{E',L') 



(4.1) 
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with the kernel 



4.2 Simplified equation for describing slow modes 



27r 27r 

= j dw\ j dw'iJ^i [r {E, L;wi),r'{E' , L';w'i)] x 



xexp!^i\^(^l[w[+l2dSi/dl'2^ - (liWi+hdSi/dh)]]. 

(4.2) 

This is a two-dimension set of integral equations relative 
to unknown functions Xh.hi-^' -^)' which are related to the 
radial part of perturbed potential by 

27r 


(4.3) 

In Eq. 1)4. we denote 

_1 (/ + fc)! {l-k)\ 

221 



(i(/_fc))!(i(/ + fc))! 



\l — k\ even, 



|;-fcl odd. 

(4.4) 



In Eq. (112 



^i{r,r') = = min(r,r'), r> = max(r, r'), 

while the function of radial action Si in (|4.2|l and (|4.3p is 



51 



dr'^J2E (/) - 2$o(r') - {h + \h\y/r'^ , 



where I = {Ii, I2, 13) and the actions I2 and I3 are related to 
the integrals of motion L and iz by L = /2 + I/3I, Lz = Is- 
The dependence r (_E, L; wi) is determined from 



r 



dr' 



y^2£- 2$o(r') - LVr'^ 



We should also keep in mind that indices li and ^2 corre- 
spond to the spatial dependence of the perturbed potential 
expanded over harmonics of the angular variables wi and 
W2, conjugate to the action variables 

^ ^ ' If""" '2E -2^o{r) - L'^/r^ dr, 



h = — (b pr dr 

2-K TV 



respectively: 

^{I,wi,W2) = (2^)-2^$,„i,(/)e'('i'"i+'='"^). 

In the case of m = 0, the dependence on the angular variable 
is absent. Let us also give an alternative form for J-i{r, r') 
(sometimes it proves to be more convenient): 

oc 

T,iry) = i2l + l)fdk-Il±10^-I^±10^, (4.5) 
J Vkr Vkr' 



where Jvjx) is a Bessel function (see 

IPolvachenko fc Shukhmanlll982l '). 



Let us assume that a massive nucleus or a black hole (with 
mass Mc) is placed at the galactic center. Moreover, we 
assume that the central mass dominates the unperturbed 
potential 'l>o(»^), so that $o(»") = $c(f) + ^G{r), where 
^c{r) = —GMc/r, $G(r) is the potential created by the 
spherical subsystem of galaxy. Then the force acting on a 
star from the central mass significantly exceeds the force 
from stars in the galactic spherical component: |$^r)| ^ 
|$Q(r)|. In this case, stellar orbits are predominately gov- 
erned by the potential of the central massive point. This 
means that we are dealing with l:l-type orbits. Due to a 
small additional potential <I>g('"), the Keplerian ellipses pre- 
cess, the precession velocity is determined by the small dif- 
ference npr{E,L) = Q,2{E,L)-ni{E,L), fJpr < r2i f72. 

An explicit expression for the precession velocity in 
terms of the potential ^cir) was found in the Section 2 
(see formula (2.8)). We shall be interested in slow modes, 
i.e., modes with frequencies of the order of precession ve- 
locities, u — 0{^lpi). Then from all items with denomina- 
tors [lu - I'l ni{E', L') - I2 Q.2{E', L')] , we keep only those 
with I'l = —I2. These denominators can be transformed into 
uj-l'iQ.i{E' ,L')-l'2^2{E' ,L') = uj -l2Vipr{E' ,L'). The 
frequency constructions in numerators of these contributions 
are equal to 

[ h ni (£', L')+l2 n2 (E', L')] OF/dE' + dF/dL' 

= I2 [flpriE', L') dP/dE' + dF/dL'] . (4.6) 

The express ion in s quare brackets is a so-called Lynden-Bell 
derivative (Lynden-BcU 1979) of the DF: 



(dF/dE)^ ripr + (dF/dL) = [dF/dL) ^ 



(4.7) 



Recall that it is defined as the derivative with respect to the 
absolute value of angular momentum, L = /2 + |/3|, with the 
adiabatic invariant, J/ = /i -f /2 + l^al , and the projection 
of angular momentum, Lz, being constant. 
Denoting 

x(-E',L)ii=-!2,!2 = (t>i^{E,L), 



n 



{E,L-E',L') = Pi^v{E,L-E',L' 



and keeping only items with I'l = — in the sum (14. ip over 
Zj, we find a "slow" equation for quantities (f)n{E,L): 



n- \^ D"' ( dE'L'dL 

n' — — / 

X P„.„,{E,L;E',L 



' n'idF/dL%, ^ 



'-n' npriE',L') 



(4.8) 



with the kernel 

P„^„,{E,L-E',L') = 

2tt 2tt 



J dwi J dw'i [r (E, L-jWi),/ {E' , L';w'i)] x 

asi/a/2-tfi)]}. 







X exp I i [71' (dSi/dl2 



2-Wi - n 



The kernel „/ can be rewritten in an explicitly real form 
if one changes the limits of integration over wi and w'l from 
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[0, 27r] to [— 7r,7r]. Then it becomes evident that r{w\, E, L) 
is the symmetric function of w\, and 

P„,„,{E,L;E',L') = 

4 /dioi / dw[:Fi [r{E,L;wi),r'{E',L';w[)] 





X cos 



(^dSi/dl2 — w'l^j cos (^dSi/dh — wi 
The quantity (f)n{E, L) can also be written in a simpler 



form 



<t>„{E,L) = 2 1 cos [n (951/9/2 x[r{E,L,w^)\ dwi. 



4.3 Simplified equation for the case of nearly 
radial orbits 

In the case of orbits with low angular momenta L we are 
interested in the set of equations (|4.8|l for slow modes al- 
lows further significant simplifications. The requirements 
imposed on the width, 5L, of localization domain of the DF 
in angular momentum were discussed in Sec. 3. 

Based on the fact that the kernel Pn„i {E, L; E' , L'), 
for nearly radial orbits, depends only on energies E and 
E' (accurate to terms quadratic in L and -L'j3 and ac- 
quires the form P„„, {E, L; E' , L') ^ P„„, (E,0; E' ,0') = 
{-!)"+"' n{E,E'), where 



n(S,S')=4 / dwi I dw[Ti[r{E,wi),r'{E',w[)] = 



(4.9) 

ni{E) = (2|g|)3/V(G M,), b{E ) = rn..4E) = GA4c/\E\, 
Vr{E,r) = V^\E\ V ~ r)/r. The unknown function 
(j)n{E,L) also depends only on E (accurate to 0{L^)): 
<prL{E,L) ^ (-!)"$(£;). Moreover, according to (g^), 
the Lynden-Bell derivative in L coincides (accurate to 
0[{Mg/Mc) L'^]) with the derivative in L, with the energy 
E constant. As a result the equation (|4.8p transforms into 
the one-dimension integral equation 





HE) 



HE) 



4nG 
21 + 1 



dE'L'dL' , 

^ u - s'Dj{E')L' ^ ' 



4.4 Model case of monoenergetic distribution 

Let us consider again the model distribution F{E,L) = 
/(!/) S{E — Eo), as in the disk case. Integrating (|4.10|) over 



E' and putting E = Eo in the resulting equation, we obtain 
the characteristic equation in the form 



21 + 1 ni{Eo) J ' Lo-s 



df{L)/dL 



vj{Eq) L 



Recall (see Sec. 2) that for the case of near-Keplerian orbits 
(i.e., orbits of the 1:1 type), the orbit precession is retrograde 
for arbitrary distributions of the potential (^cir), i.e., zu < 
0. 

For convenience, we can turn from the variable L to 
the variable v = \Qpr\ = \zu{Eo)\L = —vjL > 0. Denoting 
f{L) — f(^i^/\zu\) = fo{T^), we write 



47vG U{Eo,Eo) f , ^ 
'2l + l\MEo)mEo) J '^'^'^1.^'^^ 



dfo{v)/dv 
(4.11) 



The equation (14.111) coincides with the Eq. (2) of 
I Polvachenkol (| 1 99 1 al ) . derived immediately in the spoke ap- 
proximatiorjj The rather cumbersome derivation above pro- 
vides the basis for the spoke approach (together with the 
one for disks, in the previous Section). 

To conclude this subsection, let us note that the 
monoenergetic model under consideration corresponds to 
the specific density distribution of spherical cluster po{r) 
(and the potential 4'G(r)) and has the finite radius R — 
GMc/\Eo\. For this distribution, the quantity z^{Eo) can be 
explicitly calculated: 



. Mg 8 1 



(4.12) 



Here Mg is the total mass of the spherical cluster (recall that 
it is assumed that Mc ^ Mc). The kernel H{Eo,Eo) can 
also be calculated in the explicit form. Using the relations 
(14.51) and (14.91). we obtain 



n(£;o,£;o) 



8tv^{21 + 1) 



R 



dz 



Ci= - [Jii+^)/2{z)Ji/2{z)Y . (4.13) 



The first seven coefficients Ci calculated from (|4.13p are pre- 
sented in Table 1. After substituting these coefficients into 
(|4.11|) . we find 



1 = -Ai j vdv ^ sD'i 



dfo{u)/du 



l&n^GCi ( R 



jtuj \2GM, 
This equation can also be written in the form 



1/2 



1 = -2AX]^'A' / 

a =1 



dv 



u^d fo{u)/di' 



(4.14) 



(4.15) 



The DF fo{v) is normalized by the condition that the total 



^ The proof of this statement (and other statements concerning 

analytical properties of the functions involved, near L = 0) is ® Excluding; th e uness ential factor Qi/ir lost in r.h.s. of Eq. (2) 

omitted here. of I Polvachenkol lll991al ). 
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I 1 


2 


3 


4 


5 


6 


7 


Ci 0.135 


0.063 


0.037 


0.025 


0.018 


0.014 


0.011 



Table 1. The values of coefficients Cj. 

mass of spherical cluster is equal to Mg, i.e., J FdT = Mg- 
This gives 

j fo{u)vdv = i (27r)-^ [ti7(So)]''f^i(So)AfG. (4.16) 



4.5 Stability of the I = 1 and Z = 2 modes 

We begin with studying the stability of modes I — 1 and I = 
2. As we shall show, these modes are stable. For definiteness, 
let us take the mode I — 2. From the considerations below, 
it will be immediately obvious that they are valid for the 
mode i = 1 as well. 

We start with the equation (|4.15|l for the mode 1 = 2. 
An important point is that 1)4. 15|) for this mode (as with 
m = 1) includes only one item in the sum over s (with s = 2) 
(correspondingly, two items s = ±2 in (|4.14|) '). This follows 
from definition (14. 4p of the quantity D" and the equation 
(14451) . For 1 = 2, D'i = l, so that 

Q = -4 / 2 , 2 Q = I > 0- 4.17 



Let the DF has a beam-like form due to deficiency of stars 
with low angular momenta or, which is the same, with small 
precession velocities u. We assume that the distribution has 
only one maximum, located at = ii on the semi-axis 
^ < oo. If there is a neutral mode, at some value of 
Q = Qc, the corresponding resonance must coincide with 
the maximum of /o(j')If| i.e., uj'^ = 4u^. This means that for 
Q = Qc 

f v^\dUv)/du] 
Qc = - 2 2 du. 4.18 

J — 



Evidently, for any one-hump distribution, the right side of 
(|4.18|) is negative, as the integrand is free of singularities 
and positive everywhere. Since Q > 0, we conclude that the 
neutral mode is impossible. The absence of neutral mode 
means that the marginal value of Q, Q = Qc, which sepa- 
rates stable and unstable distributions, is also absent: each 
distribution is either stable for all values of Q or unstable 
everywhere. Since stable one-hump distributions obviously 
exist, we conclude that the mode I = 2 is always stable. 

The above considerations make it also clear that the 
conclusion about the stability of the mode / = 2 is valid only 
for the case of retrograde precession, when the quantity Q in 
Eq. (|4.17|l is positive. In the case of prograde precession, Q < 
0, so that the neutral mode (as well as the instability) exists. 
This is the well-known radial orbit instability. True, here it 
develops in a non-monotonic distribution with an empty loss 

^ The absence of the neutral mode with resonance at = is 
established trivially. 



cone, instead of the usual distributions when most stars are 
concentrated at near-radial orbits. 

However, the conclusion that the instability is utterly 
impossible in the case of retrograde precession would be pre- 
mature. The matter is that we obtain the above result con- 
cerning the mode / = 2 due to a formal reason: for this mode 
there is only one summand in the sum over s. So, indeed, the 
instability is absent for I = 2 (and I = 1 as well). However, 
for modes with 1^3, when there are at least two summands 
in that sum, the instability becomes possible under suitable 
conditions. In the following subsection, we study the mode 
/ = 3 in detail and demonstrate that the instability can 
occur here. 



4.6 The mode I = 3 

Considering the case of retrograde precession and restricting 
ourselves only to the mode I = 2, we showed (Sec. 4.5) that 
neutral modes (and consequently the instability) are absent 
for one-hump distributions. However, this is valid only in the 
special case that the mode has one resonance, as with I = 2 
OT I = 1. Recall that the proof is based on the fact that 
the resonance must then be located at the DF maximum. 
In such a situation, the characteristic relation cannot be 
satisfied as the signs of right and left sides of (|4.18|l are 
necessarily opposite, when uj equals the frequency of neutral 
mode. 

This proof, however, fails if a neutral mode has two 
(or more) resonances. Indeed, the resonances can then be 
located so that the resulting growth at one group of reso- 
nances is totally cancelled by an equal damping at another 
group. In these conditions, a neutral mode can exist. This 
means that the former group of resonances must be located 
right of the maximum while the latter group to the left. 
In the simplest case when there are only two resonances, 
we must conclude that these resonances necessarily lie on 
different sides of the maximum. Then it is hard to make 
a certain conclusion about the sign of the integrand in the 
characteristic relation that involves the principal value inte- 
grals, with a singularity at each resonance. One may hope 
therefore that we can find neutral modes (and consequently 
the instability) for / 3 when there is at least a couple of 
resonances. Now we study the possibility of neutral mode in 
the simplest suitable case of Z = 3. 

For l = 3,we have Dl = j^, = The character- 

istic equation (|4.15|l gives 

(4.19) 

Let us suggest that the neutral mode with the frequency 
Lo = luq occurs at some value of Q = Qc- For definiteness, 
we assume that For this frequency, there are two 

resonances: 

u = ui = uJo, V = V2 = \uJo. (4.20) 

^ It is apparent that with a given neutral mode of the frequency 
ujQ, a neutral mode of the frequency —loq also exists (with the 
same Q = Qc)- Thus we can seek a frequency a;o squared for 
neutral mode. 
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Obviously, the resonance corresponding to smaller f (i.e., 
V = uJo/3), must lie to the left of the maximum of the func- 
tion foii^) (we denote its position it), while the resonance 
corresponding to larger v (i.e., v — uio), must lie to the right 
of the maximum: ^ luq < u < uiq. Bypassing the singular- 
ity in the complex plane from below and equating the 
imaginary part of the full integral (|4.19p to zero, we find 

/o(^o) + |./-o(|^o) =0. (4.21) 

Here we should explain that direction of bypassing in the 
complex plane coincides with that in the complex plane 
1/ (i.e., it is from below) because of uo > 0. The condition 
1)4. 21|) expresses the balance between growth at one of reso- 
nances and damping at the other. Eq. (|4.2ip determines the 
frequency of the neutral mode which exists when Q — Qc, 
the latter found from the condition 

Eq. (|4.22|) involves the principal value integrals. The 
pair of equations, (I4.2ip and (I4.22p . determines the fre- 
quency of the neutral mode and the critical value of param- 
eter Q, Q = Qc such that the system has a neutral mode. 
This is in fact the condition on some DF parameter (say, 
dispersion of precession velocities, ut), or on mass, AIq, of 
spherical component (i.e., on value of its self-gravitation). 
Under this condition, the system is at the stability bound- 
ary. If the parameter deviates from its critical value in a 
certain direction, the system becomes unstable. This direc- 
tion is yet to be determined. 

It is not evident beforehand that the right side of Eq. 
(I4.22P will be positive for potential neutral modes. (Their 
frequencies are determined from Eq. (I4.2ip . It is easy to un- 
derstand that for one-hump distributions, a suitable solution 
of this equation always exists.) 

Moreover, we showed above for the mode I — 2 (in fact, 
for the mode Z = 1 as well) that in principle the neutral 
modes are absent as the resulting value of Q turns out to be 
negative. 

In order to clarify the possibility of neutral modes 
Z 3 in more detail, we consider the series of specific 
models in the form of one-hump distributions /q"' {ly) = 
Nn{v'^)^ exp (j-v^ /v't) , where A'^„ is the normalized coeffi- 
cient, n = 1,2,... . On the assumption that the DF is nor- 
malized to some total mass M by the condition 

roc 

/ uduft^=M, (4.23) 
Jo 

then Nn = [2M/(!/|.)"+^] (n!)"^ Dimensionless variables 
are introduced for convenience: x — jv^, — ujq/v^. 
Besides, we shall use the function fn(x) — /g"' (ut^/x) = 
Nnivx)" x" , instead of the function fi"\i')- In new 
designations, the set of equations (|4.2ip and (|4.22p takes 
the form 

fUxo) + 4f.fL{^xo)^0, (4.24) 

oo 

Q^^_Idxx^M-l(^ + ^), (4.25) 
J ax \xo — X Xa — yx J 



where 

fUx) = (iV„ 4") a^""'(n - x) e-\ (4.26) 



71 








1 


1.66 


3.44 


8.85 


2 


2.27 


5.343339 


17.99 


3 


3.09 


7.733736 


26.99 


4 


4.03 


10.17 


35.99 


5 


5.01 


12.63 


44.99 


6 




15.09 




7 




17.55 





Table 2. The values of squared dimensionless frequencies — can- 
didates to neutral mode. 

One can see from (|4.26p that the maximum of the function 
fn{x) is located at x = n, so that (|4.20p gives the condition 
on the value of dimensionless frequency squared: n < xq < 
9n. The value xq is found from Eq. (|4.24p . that takes the 
form of a transcendental equation 

(9n - a;o) - i 3^"+^ exp (-§ lo) {xo - n) = 0. (4.27) 

It is evident that at least one root always exists as the left 
side of Eq. (|4.27p has opposite signs at the ends of the in- 
terval under consideration. More detailed calculations (or 
plotting the left side) show that there are actually three 
roots satisfying the condition n < xq < 9n, for each n. All 
roots are candidates for a possible neutral mode. The nu- 
merical solution of Eq. (|4.27p . for n = 1, 7, gives squared 
dimensionless frequencies of neutral modes listed in Table 2. 

We next substitute the obtained values of xo into the 
relation (|4.25p that determines Qc- Let us represent it in the 
form convenient for numerical calculations. To do this, we 
introduce the function gn{x): 

g„(x) — e [{x — t) (X — t — n)\, x > 0. 

— oo 

Using this function, Eq. (I4.25P can be written as 

Qc'^ = ^^(a^o) + I gn i^xo) , 

where we denoted Q'"' = Q/K" iV„] = n! {Q/2M) v^, and 
the quantity M (see (|4.23p l is determined by the normaliza- 
tion requirement for the monoenergetic model, (|4.16p . i.e., 

M = |(27r)-' ^x7(£;o)'f^i (£())Mg. 

Substituting all three roots x'^^ , Xq^' , x^^^ into (|4.25p . we ob- 
tain the following results. 

(i) . For the DF with n = 1, the right side of (|T25ll 
turns out to be negative for each potential neutral mode. 
This means that the mode I = 3 has no neutral modes (con- 
sequently, the mode is stable) for the model n = 1. 

(ii) . For the models with n ^ 2, there is one neutral 
mode. It corresponds to the middle root xq (in Table 2, it 
is the root x^^^). Only for this root, the right side of (|4.25p 
turns out to be positive. This means that these models can 
be unstable if the parameter Q differs from the critical value 
we found. Recall that the deviation direction of Q is yet to be 
determined. The critical values Qc"' are presented in Table 
2. 

Fig.|S]shows the resonances for the neutral mode in the 
model with n = 2, for illustration. As is seen in Table 1, the 



© 2006 RAS, MNRAS 000.[Tlf26l 



16 E. V. Polyachenko et al. 




Figure 8. Positions of resonances for the mode Z = 3 for DF 
/„_2(!^) oc u"^ exp{— i^^/i^y). The frequency of a neutral mode 
is ujQ = Vs. 34 ut = 2.31ut- The DF is maximum in the point 

U = U = \f^VT- 



dimensionless frequency ujq/vt is equal to V5.34 ~ 2.31. So 
the resonances lie at v/vt ~ u-io/(?)1't) ~ 0.77 and u/ut = 
wo/j^T = 2.31. 

We are thus convinced that the neutral modes I 3 
(and consequently the loss-cone instability) indeed exist for 
one-hump distributions for suitable parameters (n ^ 2). 



4.7 The perturbation theory near Q — Q 
instability criterion. 



(n) 



The 



Let us impose an increment SQ on the parameter Q and cal- 
culate a correction to the squared dimensionless frequency, 
5xo, using the perturbation theory. Then the instability 
corresponds to the positive imaginary part of Sxo- Indeed, 
Sxo = S{lo'^)/i'^ = 2uJoSio/u^, and the signs of imaginary 
parts of Sxq and Suj coincide as luq > by agreement. 

Let us write the characteristic equation in the form 
(|4.19|) . using the new variables x = v^/u^, xo = luq/v^ 
and a new function fn{x) — x"e~^. 

oo _ 

dfnix) ( 1 



dx X ■ 



dx 



+ 



5/3 



X — Xo 



^XO 



(4.28) 



where /„ (x) = x" ^{n - 



We find 



SQ 



(n) 



5x0 I dx X 





df ujx) 
dx 



+ 



5/27 



{x-xoy (^r^^^xof 



where the integration involves bypassing the singularities 
from below. Now we apply the relation useful for integrating 
the expressions with peculiarities of the type {x—xo)~ (with 
indentation due to the singularity in the complex plane): 

b b 

F(x) , f F(x) 

^ ^ ■dx = FP ' ^ ' 



{x — xo)' 



{x — xo)' 



■dx + i-RF'{xo), (4.29) 



where a < xo < b and FP means "Finite Part" of the inte- 
gral. It is easy to reproduce the derivation of this formula 



when we recall the definition: 



FP 



f{x)dx 
{x - Xo)'^ 



Yim\ 

e-.0L 



+ 



f{x)dx 2f{xo) 



fix)dx 

{x - Xo)'^ ' J {x - Xo)'^ € 

XQ + E 



where a < xo < b, and the function f{x) is regular at xo- 
Evidently, the direction of bypassing (either from below or 
above) has no influence on the real part of the result (i.e., 
the value of the FP integral). Changing the direction of in- 
dentation, we change only the sign of the imaginary part in 
(|4.30l) . [Note that the concept of FP integral is well-known 
in hydrodynamics. It is widely used in problems related to 
the so-called critical layer, i.e., a narrow domain near the 
resonance of the wave and the shear flow of fluid (see, e.g., 
Hickernell 1984, or Shukhman 1991).] We obtain 



rf/n(a 



+ FP 



dx X 



dx 
dx 



+ 



[{x-. 



+ 



5/27 



[x-xo)'^ ( 



: xo) 



(4.30) 



The imaginary part (|4.29p can be simplified using the re- 
lation (|4.24p . reflecting the balance between growth and 
damping for a neutral mode. Then we obtain 



(JQ^"-* = 5xo\ i-KXo [f^lixo) + 2I3 fn {k ^0) 



+ FP / dx 
Jo 



xfU^)\-, — - 

LIx — X 



.{x - Xo)'^ 
Writing (|i3T1l in the form 5Q^"^ 



+ 



5/27 



(4.31) 



[x-^xoY 
(A„ + iB„) Sxo, where 



00 



A„ ^FP I dx X /„ (x) 




+ 



5/27 



(x-xo)^ ( 



X - -^xo 



Tvxo (a;o) +233/" (| 2;o)] 



we find for real and imaginary parts of Sxo'- 



Re (Sxo) — 



An 



SQ^"\ Im{Sxo) 



Bn 



A'^ + B?. 



5Q 
(4.32) 



We see that the stability criterion is determined only by the 
sign of Bn- Calculating the growth rate itself requires the 
values of both quantities. An and Bn- Using the functions 

Un(x) = X fn"{x) — x"~^ — 1 — x) {ti — x) — x^ , 

X 

fe„(a:)=e-"FP j ^ {x - t)''{n - x + t), x > 0, 

— oo 

the result can be written in a compact form 

An = hn{xo) + ^ hn (| Xo) , Bn ^ TV ^Un{xo) + §f U„ a^o)] ■ 

The values of An and Bn are calculated numerically, for the 
values of xo found above. They are presented in Table 3. It 
is seen from this Table that Bn is positive for all values of 
n. Consequently, the instability occurs when (5Q("' < 0, or 
Q'"' < Qc"\ where the critical values of Qc"' are presented 



© 2006 RAS, MNRAS 000.[Tlf26l 



Gravitational Loss- Cone Instability in Stellar Systems with Retrograde Orbit Precession 17 



n 


xo 










2 


5.343339 


0.542582 


0.2713 


0.188022 


0.732157 


3 


7.733736 


3.486635 


0.5811 


0.218900 


1.884489 


4 


10.17 


15.6856 


0.6536 


0.08 


5.46 


5 


12.63 


74.7372 


0.6228 


-1.40 


18.32 


6 


15.09 


399.7582 


0.5552 


-11.86 


70.95 


7 


17.55 


2425.3916 


0.4812 


-83.65 


312.89 



Table 3. Dimcnsionless frequency squared xg, the value of critical 
parameter Qc'^ and the values of A„ and B„ in the expression 
for complex frequency squared Sxq = S{ui'^ /u^). The unperturbed 

'■'")/n(x), /„(X) 



DF is U{x) = {Nr,yl 
there are no neutral modes.) 



' exp {—x). (For n = 1 



in Table 3. If we recall the definition of Q<"' = n\{Q/2M)vl, 
the instability condition of the mode I — 3 for the monoen- 
ergetic model can be reformulated in a form of criterion on 
the ratio of the angular momentum dispersion, Lt ~ vt jvj, 
to the angular momentum of a star with the same energy 

Eq in a circular orbit, 



\ GMcR: 



Lt TV J 6C3 I 



C3 « 0.0373 



(4.33) 



Particularly, for the model with n = 2, when Qc ~ 
0.5426, we obtain from Lr/icirc < 0.193, and for the 

model with n — when Q^''' ~ 3.4866 we obtain Lt/ Ldrc < 
0.283. Note that the criterion in such a form does not involve 
mass of the spherical component, Mg0 

When the supercriticality of SQ is not small, the pertur- 
bation theory above does not allow us to calculate the com- 
plex eigenfrequency. In this case, the characteristic equation 
(|4.28|) was solved numerically, for values of n = 2, ... ,5. The 
qualitative behavior of real and imaginary parts of the fre- 
quency is similar for all calculated cases. So we restrict our 
illustrations only to the model with n — 3 (see Fig.|9]). Note 
also that the results of computations, for small deviations 
from the stability boundary, coincide with the asymptotic 
results obtained using the perturbation theory (|4.32|) . 
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Figure 9. Loss-cone instability in spherical system, a — DF 
over precession angular velocity /n=3 = x^ exp(— x). The dash 
lines indicate positions of resonances for neutral mode. (Squared 
dimensionless frequency is xq = 7.733736.); 6 — Behavior of 
squared dimensionless frequency xq of unstable mode. The values 
of Re(xo) and Im(a;o) are shown against parameter Q. 



form 



First we represent our characteristic equation in the 



4.8 General instability criterion and study of 
specific distributions 

The above results, based on the neutral mode approach, 
can also be obtained using a suitable analogue by means of 
the well-known Penrose - Nyquist theorem (Penrose 1960, 
Mikhailovsky 1974). Recall that this theorem is widely used 
in plasma physics. Employing the theorem helped to estab- 
lish numerous general results in the theory of plasma insta- 
bilities. 



* Since these critical values (-I/T/£'circ)crit proved to be not too 
small, it means that in accepted "spoke approximation" (where 
it is supposed that L/L^irc <S 1) instability criterion is know- 
ingly satisfied. However the rigorous calculation of the instability 
boundary in terms of parameter Lt / icirc requires an exit from a 
framework of this approximation. This problem is under consid- 
eration now. 



= Ci / dxx 



.df{x) 



dx 



E 



Df 



z/ s 



(4.34) 



where the quantity Q is independent of I. In the case of 
retrograde precession we are interested in, Q > 0. Recall 
that Sinin is cqual to 1 or 2 depending on evenness oil, x = 
is the precession angular velocity squared, in dimensionless 
units (the units of vt, where vt is the precession velocity 
dispersion, is common), f{x) is the unperturbed DF, z = 
is the frequency squared, in the same dimensionless units. 

We are also reminded that the singularity of the integral 
in the right side of (|4.34p . for z on the real axis, must be 
bypassed from below if Re (uj) > 0, and above if Re {u) < 0. 
From this indentation rule and the form of Eq. (|4.34p . it 
immediately follows that complex unstable roots, zo, form 
pairs: if z = zo — a + ib, a ^ 0, b > - the root of 
Eq. (|4.34|l . and the corresponding eigenfrequency is cij = 
too — a + i P, {a > 0, /3 > 0), the complex conjugate 
root, z = Zo = a ~ ib, also satisfies Eq. (I4.34|) and describes 
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the mode with the same growth rate, but opposite sign of 
frequency, to — —uo = —a + 1/3. 

Note that the form of Eq. (|4.34|) shows that the ape- 
riodic instability, for which Re (20) < 0, lm{zo) = 0, i.e., 
a = 0, is absent here. Indeed, putting zo = —\zo\ < in 
(|4.34|) and integrating by parts, it is easy to verify that the 
right side of (|4.34|l is negative. By the way, this is a fur- 
ther distinction of Eq. (|4.34p from the disk characteristic 
equation that aUows aperiodic instabilities for symmetric 
distributions /(i'). 

Eq. (|4.34|l involving [i (I + 1)] items can be reduced to 
a one-term equation. Indeed, the substitution 

I 

F^'\x) = Ci Dtfix/s"). (4.35) 

S — Smin 

transforms (14.34p into the equation 

Q^Jd.J-I^^Ml±. (4.36) 
J X- z 



Thus, for arbitrary values of Z, we obtain a one-item equa- 
tion (|4.36p similar to that for the mode I = 2 (or I — 1). 
However, now the integral involves the function F'^'^x), in- 
stead of the initial function f{x) (with one maximum and 
tending to zero at a; = and infinity). Starting with modes 
I = 3, the new function, F'^''>{x), can have minima. It is easy 
to understand that the frequencies (candidates for a neu- 
tral mode) calculated in subsection 4.6 from the condition 
of balance between growth and damping at resonances on 
different sides of the maximum of the initial function f{x), 
are precisely those coinciding with the extrema of the new 
function, F^'\x), i.e., z = Xj, F'ixj) — 0. Correlating this 
with earlier results for disks, and also with the sphere DF 
fn{x) = x" exp (—a;), we become convinced that the largest 
(more often, the only) positive critical value of Qc for the 
neutral mode should necessarily be related to a minimum 
of the DF F^'\x). As an illustration of this statement. Fig. 
10 shows three functions F^ ^ (x) of this series, for the mode 
I — 3. From the figure (and also Table 2), it is seen that 
only the central of these three extrema, i.e., the minimum, 
gives rise to the neutral mode with positive values of Qc (for 
n ^ 1). 

Thus, we see that the availability of minima is of funda- 
mental importance for the existence of neutral modes with 
positive Qc (and consequently for instability). We have al- 
ready seen that for / = 1 and I = 2, when the DF F''^'>{x) 
coincides with the initial DF f{x), so that the former has no 
minima, the corresponding neutral modes (and instability) 
are absent. 

Recall that the possible-in-principle neutral mode with 
z — corresponding to the resonance at the minimum a; = 0, 
has Qc ~ 0, so that it cannot be assumed to be a candidate 
for a neutral mode with the property required for instability 
(Qc > 0). Therein lies a fundamental difi'erence from disks 
where any two-hump distribution with a zero minimum at 
1/ = always has the neutral mode with Qc > at the 
minimum. So such a distribution is always unstable (when 
Q < Qc) independently of other DF details. 

We have a right to expect a neutral mode (and insta- 
bility) related to the minimum, for I J5 3 only. In fact, it has 
already been demonstrated above using a somewhat more 




5 10 15 20 25 30 35 40 



Figure 10. — Function Fn\x) for a model f(x) = x"e ^ with 
i = 3 (n = 1, 2 and 3). 

cumbersome method. Besides, a question remains unsolved 
in the approach we apply: why are not all distributions 
with one maximum that vanishes at the ends of the posi- 
tive semi-axis ^ a; < 00 generally unstable, even though 
1^3. Empirically, by considering various series of distribu- 
tion functions, we found the qualitative instability condition. 
Its rough formulation is: the instability is possible if the DF 
function is well-localized around its maximum. 

Now a possibility appears for a more rigorous (in 
fact, quantitative) formulation of the instability condition. 
Though Eq. (|4.36|l differs from the equation 

oc 

Q= f (4.37) 

J V - c 
— 00 

(where c — co/k is the complex phase velocity, Q — /ljo > 
0, k is the wave number, ujq — i-rmoe ^ /m is the plasma 
frequency squared) , for which .Penros^ (|l960l ) obtained his 
well-known criterion, here we also can obtain an analogous 
criterion ~ i.e. a counterpart of the Penrose - Nyquist crite- 
rion, for our equation (I4.36|l . First we formulate it in terms 
of neutral modes. 

Theorem. The distribution F'^''\x) is stable if neutral 
modes corresponding to minima of F^'\x) are absent. Al- 
ternatively, if at least one neutral mode corresponding to a 
minimum occurs, then a sufficiently small Q always exists, 
for which the system will be unstable relative to perturba- 
tions with a given I. 

The instability condition for any / follows immedi- 
ately from the theorem. Indeed, if for at least one of I 
{I — 1,2, . . .), a neutral mode exists for the corresponding 
distribution F'^^\x), then such a sufficiently small Q occurs, 
for which the system is unstable. 

It is useful to give another formulation of t he theo- 
rem with a maximally possible similarity to that of |Penros3 
(1960) for Eq. (1133. 

Theorem (another formulation). The distribution 



© 2006 RAS, MNRAS 000.[Tlf26l 



Gravitational Loss- Cone Instability in Stellar Systems with Retrograde Orbit Precession 19 



F^'\x) is stable if and only if for all points Xj, at which 
the modified DF F^''>{x) has a minimum (i.e., F'{xj) = 
0, F"{xj > 0), the condition 

f, dF<-'\x)/dx „ 
Qc= / dxx ^— < 0. (4.38) 

J X- Xj 



is met. Conversely, if at least for one minimum the opposite 
inequality is satisfied, then such a sufficiently small Q exists, 
for which the system is unstable for perturbations with a 
given I. 

The proof of this theorem can be found in Appendix 
B, but now we discuss a correlation between the instability 
condition following from this criterion and our qualitative 
condition formulated above. 

From the results obtained for disks, we know that under 
sufficiently deep minimum, the corresponding Qc can be- 
come positive. Thus, it can rigorously be shown that Qc > 
(with a finite margin) in the limit when the minimum is 
exactly equal to zero. 

Indeed, let us assume that Xj is the position of mini- 
mum, at which F{xj) = 0, F'{xj) = 0. Then we obtain by 
integrating in (|4.38|l by parts 

OO OQ 

f , dF'^^Ux)/dx f , dF'^^\x)/dx 
Qc= I dxx = Xj I dx = 

J X Xj J X Xj 



7 , F^'Hx) 
= Xj dx ^—^ > 0. 



It becomes impossible to prove in a similar manner that 
the integral is also positive for a non-zero minimum. Actu- 
ally, it can have any sign. However, positive contributions 
into the integral increase the closer the minimum is to zero. 
So the integral should eventually become positive with in- 
creasing depth of minimum. 

In light of this fact, it becomes clear that our qualita- 
tive instability condition means that a minimum of F''-'(x) 
is sufRciently deep, so that the related neutral mode occurs. 
This becomes evident when we consider how the function 
^(''(a;) is buih from the initial DF f{x). The instability 
is unavailable if a minimum is absent or is not sufficiently 
deep. In Fig. 10, the functions Fn\x) constructed from the 
functions fn(x) = e~^, for the mode I — 3, are for conve- 
nience calibrated so that the value of the highest maximum 
is equal to unity for all values of n. We see that the only min- 
imum becomes deeper with increasing n. This is in complete 
agreement with the results of Sec. 4.7 where we found that 
the mode i = 3 is stable for n — 1 and unstable for n ^ 2. 
(For completeness we additionally checked that instability 
emerges when n > 1.55). Moreover, we checked that the 
instability in the model with n = 1 is also absent for any I. 

The formulated criterion allows a purposeful search for 
such DF f{x) which gives a new function F{x) with a mini- 
mum capable of 'generating" a neutral mode. In other words, 
the integral (f^35)l , 

OC OO 

^ f, dF^'\x)/dx f, FW(x)-FW(a,) 
Qc = dxx = Xj / dx ,„ , 

J X-Xj J [x-XjY 



(4.39) 

is positive, so that the instability occurs when Q < Qc- It 



turns out that suitable distrib utions are known in plasma 
physics. In particular, Penros^ l|l960l ) has pointed to a case 
of such a distribution. Namely, he demonstrated that a 
plasma distribution becomes unstable provided this distri- 
bution has a sufficiently sharp minimum (then the Penrose 
integral similar to (|4.39[) becomes positive). For instance, 
such a minimum appears at the electron DF when a suffi- 
ciently cool electron beam is injected into the Maxwellian 
plasma, provided the beam velocity is larger than the elec- 
tron thermal velocity of main plasma. 

It is interesting that distributions with sharp minima 
appear in our problem quite naturally. Indeed, let us assume 
that the star distribution with respect to angular momenta 
(or, which is the same, to precession angular velocities) is 
Gaussian. In terms of the variable x = v'^ /v^, it has the form 
f{x) = Ne~^, < X < OO, A'" is the normalized constant. 
Now we suggest that the stars enclosed by the loss cone elude 
the distribution so that the resulting distribution arises: 

f{x) = N H{x-a)e''', (4.40) 

where H{t) is the Heaviside step function. Any physically 
admissible distribution should of course be smooth. So, in- 
stead of discontinuous function (|4.40|l . we assume a nearly 
identical (but continuous and smooth) distribution 

fix) = A7^a(a;)e~", 

TZa{x) = i [tanh (^^-^) + tanh Q)j , (4.41) 

where 5 ^ a in the "cutting factor" TZa(x). The function 

F{x,a) = fix) + ^ f{\x)\ 

for the mode / = 3, corresponding to the DF (14.411) with 
a = 0.1, a = 0.01 and 5 — ^ a, is shown in Fig. 11 (a). It 
is seen that the distribution Fix) has only one minimum, 
this minimum being sufficiently sharp. Direct calculations 
show that the neutral mode (Qc > 0) corresponding to the 
minimum occurs under an arbitrarily small size of a "slot" 
(i.e., a value of a). Fig. 11 (6) shows the marginal curve on 
the plane (Q — a) , where the modes with arbitrary values of 
I are taken into account. It is seen that for not too large val- 
ues of a the boundary of instability is nevertheless actually 
determined by the ffrst unstable mode Z = 3 only. 

Thus, we see that an empty loss cone, even if it is very 
narrow, inevitably leads to the instability, for suitable distri- 
butions (i.e., the dispersion ut is less than the critical value 
(i^T)c determined by the parameter Qc.) 

5 DISCUSSION 

First we list the results of the paper. 

1. The paper presents a systematic derivation, from gen- 
eral linearized Vlasov equations (written in the action-angle 
variables), of simple characteristic equations for small per- 
turbations in disk and spherical stellar systems with near- 
radial orbits. 

2. On the one hand, our analysis of these characteristic 
equations conffrms the presenc e (alr eady discussed earlier, 
Polyachenko 1991b, Tremaine [iooi) of gravitational loss- 
cone instability in disks. On the other hand, we succeeded 
to prove, in this paper, a possibility of this instability in 
spherical clusters. 
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Figure 11. (a) - Function F(x,a) = C3 [3 /(x) + 5 /(| x)] for 
the model (8) with I = 3 and a = 0.01, a = 0.1 (5 = a/20). (6)- 
the marginal curve on the (Q — a) - plane, calculated for all I. 



3. It is shown that the physical reason of the instabil- 
ity under consideration is escape of stars through the loss 
cone due to destruction of stars with sufficiently low angu- 
lar momenta. As a result, the DF in angular momenta will 
assume an unstable ("beam-like") form. This is very similar 
to the situation in plasma traps (mirror machines) where 
plasma particles with low transversal (relative to the axis 
of the trap) velocities escape those systems. For this reason, 
distribution over these transversal velocities also becomes 
"beam-like" , so that the classical loss-cone instability devel- 
ops. 

4. We highlight retrograde precession of orbits as a nec- 
essary condition for the gravitational loss-cone instability. 
Expressions for precession velocity both in non-singular and 
near-Keplerian potentials are derived. In particular, they 
helped to obtain conditions for the the precession to become 
retrograde. 

5. While deriving the characteristic equations, we justify 
the obvious (and very convenient for practical use) rotating- 
spokes approximation. 

6. For analyzing the characteristic equations, a specific 
method is developed. It is based on preliminary search of 
neutral modes. 

7. We also developed a method based on generaliza- 
tion of the plasma (and, in fact, gravitating disk) Penrose - 
Nyquist theorem that establishes the criterion of stability 
and instability. First, from an initial DF of stars in a spher- 



ical cluster, /, we turn to another, effective DF, F. The lat- 
ter is constructed from / according to a simple recipe (see 
the beginning of Subsec. 4.8). Using this new function, the 
many-term characteristic equation describing perturbations 
in spherical systems reduces to the simplest one-term disk- 
like equation. In turn this allows us to formulate and prove 
the above-mentioned generalization of the Penrose - Nyquist 
criterion. 

8. It is shown that this criterion allows us to justify 
the following qualitative criterion of the gravitational loss- 
cone instability: the instability is possible if the DF is well- 
localized about its maximum. Using the criterion, we can 
perform a purposeful search of unstable distributions. In 
particular, we succeeded in proving an empty loss cone, even 
if very narrow, to be able to lead to instability. 

Let us remind that for spherical near-Keplerian sys- 
tems, Tremaine does not find the instability, although his 
integral equation (55) is equivalent to our "slow" equation 
(4.8), and the problem definition is very similar to ours: 
To study the stability of spherical models against compar- 
atively slow perturbations (with typical times of the order 
of inverse precession frequency) provided that DF possesses 
positive derivative dF/dL > at low angular momentum 
due to the loss-cone. 

Still, t here are tw o considerable differences. The first 
one is that iTremaind (|2005) has used Goodman's (1988) 
criterion, which is fundamentally restricted to the modes 
/ ^ 2. The second one is that Tremaine studied only mono- 
tone DF (from radial orbits L = up to circular orbits 
L = Lcirc). In that case the variational principle takes place, 
which claims that squares of eigen frequencies should be real 
(Im(ci;^) = 0). S o the unstable modes are aperiodic, if any. 
I Tremaine! (|2005l ) has shown that for the models considered 
(log-models, see (98) in his paper), stabilizing and destabi- 
lizing contributions cancel each other for / = 1 leading to 
neutrally stable lopsided mode uj = 00 while for ? = 2 the 
stabilizing contribution dominates. 

In our case non-monotone DFs violate the variational 
principle, so that the squares of eigen frequencies should 
not be real. Moreover, we have shown (Sect. 4.8) that the 
aperiodic instability is impossible, i.e. the eigenmodes are 
oscillating, Re(a;) 7^ 0. We have found the instability at i 3 
for spherical models, if DFs have maximum somewhere in 
the region < i < LcircH It is plausible that we deal with 
instabilities of somewhat different nature. 

By considering two-humped models in Sect. 3, we 
demonstrated that the gravitational loss-cone instability can 
arise in disks with nearly radial orbits. One can show that 
two-humped DFs is crucial for the instability. Such distribu- 
tions arise naturally from originally Maxwell-like DFs after 
consumption of stars with low angular momenta by a black 
hole. The instability takes plac e for arbitrary azimuthal 
number m. Meanwhile ITremaind (|2005l ) has considered dis- 



^ The lopsided zero mode cannot be found with our spoke-orbit 
approximation (4.15). To obtain the mode, one should hold terms 
of the order of O(L^) in the expansion of functions in (4.8). 

Strictly speaking, in the spoke-orbit approach we use here, DFs 
describe almost radial orbits, L <^ icirc- However, we believe that 
presence of maximum, rather than high degree of orbit elongation 
is of fundamental importance for the instability. 
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tributions in the interval — Lcirc < L < Lcirc, which is sym- 
metric around L — 0. On the example of lopsided mode 
Z = 1, he has shown that disk is generally unstable. Note 
that condition (3.11) is implicitly susg ested in hi s deriv a- 
tion of instability criterion (Eq. (114) in IXremaind (|2005l )). 

For the final conclusion on which of two differences is 
decisive for the instability, i.e. monotonity or mode number 
I ^ 3, some additional work is needed. We hope to attack 
this problem in future using our "slow" equation, derived in 
Sect. 4. The problem seems to be of importance since in some 
numerical models (e.g., Cohn and Kulsrud, 1978) collisions 
result in the establishment of distributions growing mono- 
tonically (from the loss-cone boundary up to circular orbits) . 
It is likely that such distributi ons will prove to be s table (see 
recent N-body experiments bv lBerczik et al\ l|2005l 'l'). On the 
other hand, we also emphasize the fact that such a "near- 
isotropy" observed in some numerical simulations is not an 
unambiguously established fact at present. Not to mention 
that distributions with prevalent elongated orbits can be 
quite natural in some circumstances (e.g., periods between 
their formations due to, say, collisionless collapse, and the 
moment of relaxation)^ 

In the paper by iBerczik et all l|2005l 'l N-body simula- 
tions have been used to model the loss cone in the vicinity 
of the (binary) black hole. The model allowed them to con- 
struct the lose cone which was maintained in nearly-empty 
state during simulations. However the black hole feeding 
rate was at the level consistent with standard collisionally- 
repopulated loss-cone theory. The discrepancy between their 
and our results is probably explained by diffe rence in char- 
acter of distribution over angular momentum. iBerczik et al\ 
l|2005t ) have assumed initially homogeneous system, which 
leads to the near-isotropic one in the course of collision 
evolution (more exactly, to the distribution slightly grow- 
ing with L if to ignore a small region of empty loss-cone), 
whereeis we have proved instability for systems with nearly 
radial orbits. Yet, this N-body modelling can be a weighty 
argument in favour of stability of monotone distribution 
functions. 

In conclusion we present preliminary estimations of effi- 
ciency of the collective mechanism under consideration. For 
the most interesting, ne ar-Ke plerian, case, such estimates 
were made by Tremaine l|2005l '). and we use these estimates 
below. 

There are several characteristic time scales. The first is 
the dynamical time, tdyn ~ ~ {E? /GMc)^^^ , where R 
is the typical orbital radius, Mc is the central point mass . 
The orbit precession determines, using Tremaine's (120051 ) 
terminology, the secular time scale. 
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Me 
Mg 



tdy 



Nam ' 



where Mq, Nq and m is the cluster mass, the number of 
stars and the mass of one star, respectively. The gravita- 
tional loss-cone instability develops precisely on this time 
scale (cf. the formula (|4.12p for the precession velocity in 
our monoenergetic model, 7 = Im (uj) ~ fipr — voL ~ 
(Ma/Mc) {L/R^) ~ Q. {Ma/ Ma).) The next important time 
scale defines a period of collision relaxation. 



These three time scales are well-known. iTremaind (|2005D 
introduces another (less known) time scale, citing his paper 
iRauch fc Tremain|[l996l ') - the time of resonance relaxation 
of angular momenta, 

R^'^Ml'^ 



Gi/2 



^dyn 



R"H{t'^ 

G^/'^mMa 



tdy 
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niMa 



tdy 



m; 



For near-Keplerian systems (when Mg ^ Mc, Ng ^ 
1), these four time scales are highly separated: 

^dyn ^scc ^rcs ^rclax. 

Thus, according to these estimates by Tremaine, the insta- 
bility should grow faster than collisional (and resonant) re- 
laxation, whether or not a cluster mass is small. 

Note, however, that the estimates of time scales pre- 
sented here are insufficient to claim that the collective mech- 
anisms under consideration should dominate. There is a need 
to calculate and compare the star fluxes onto the black hole. 
In this connection, we should remind the reader that so far 
we only attempted to prove that the instability is possible 
in principle. 
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the components Si, S2 and S3 are 

r 

Si= f dr'^2E {I) - 2$o(r') - (J2 + Ihlf/r''^ 



S2^ J Pede'^ J sj{l2 + \h\y -Il/sin' 6' dO' , 

^3 = J h dip' = I3 if. 


For convenience, we accepted a symbolical, "vector", 
designation / = {Ii, 12,13). We have 



wi = dS/dh = dSi/dh = Qi 



W2 =dSi/dl2+ de' 



r 

I 



dr' 



y'2£;-2$o(r')-iVr'2 



I2 + 1/3 1 



\/(/2 + |/3|)2-/37sin' e' 



APPENDIX A. Derivation of the integral 
equation for perturbations in a 
spherically-symmetric gravitating system in 
terms of the action-angle formalism 

1. The action-angle variables in a 
spherically-symmetric potential 

Let us recall the action-angle variables in a spherically- 
symmetric potential $o(r). 
The action variables: 



Ji = -!- prdr = - 
An IT 



dSi 1 9/2 + arccos 



^2E-2<i>o{r)-L^/r^dr, {Al) 



/2 = — * pede = - / d6i^L2-Li/sin^6l = L-lL^ 



^3 = ^ f p^dtp = Lz 



{A2) 
(A3) 



Here E = ^Vr + \ I? jr^ -f <I>o(''') is the particle energy, L = 
r-^iig -I- = y^Pa + pj/sin^ Q is its angular momentum 
magnitude, and = r sin Ov^is a, projection of the angular 
momentum on the axis z. The angle 9q is defined as sin^ 60 = 
LI/L^, and the generalized impulses are defined as follows 
Pr — r — Vr, pe = = r vg, p^p = sin^ 9 (f — rsinOv^. 
Note that it follows from (A2) and (A3) that L = I2 + 
\h\, Lz = h. 

The angular variables. 

By definition angular variables w-i, W2, W3 are Wi = 
dS/dli. The function of action S in a spherically-symmetric 
potential is known to allow for separating the variables and 
can be written as a sum S{I; r, 9, if) — Si + S2 + S3, where 



/ cos 9 
\ cos 6*0 



W3 = ^Sl/^I3 + ^S2/^h + 'P■ 
Action variables are integrals of motion, and angular vari- 
ables linearly depend on time: Wi{t) = 'Wi{0) + ^li{I) t, where 
frequencies f2j(/) are Qj — dE[I\, I2, 13) / dij. 

2. The solution of the kinetic equation, calculation 
of perturbed density and derivation of the integral 
equation 

The perturbation of the DF /i is easily obtained from 
the kinetic equation if we write it down in terms of action- 
angle variables: 

^ = ^+n,^ = — — 
dt dt dwi dli dwi 

We have 

(27r)3 ^ 

(A4) 

This involves summation over a repeating index j = 1, 2, 3. 
In what follows the background DF F is supposed to be 
dependent on E and L only. The function ^ii, 12, 13 (-f) a-P- 
pearing in (A4) is 

2ir 2ir 2ir 

$ii,i2,i3(-f) = j dwi J dw2 J dw3 ^{I,w) exp{ilw), 

000 

(A5) 

where another symbolical "vector" designation are intro- 
duced for brevity: w = {wi,W2,W3) and I = {li, 12,13). In 
Eq. (A5) the function "!>(/, w) represents the perturbation 
of potential (without the factor e~*"') expressed in variables 
(/, w). We shall choose this function in the form 



<i>{r,9,if)^xir) Piicos9), 



{A6) 
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where Pi{x) is Legendre polynomial. In the main text we 
have already given arguments why we may confine our- 
selves to the case m = 0, without considering perturba- 
tions with a more general angular structure of the type 
Yi"'{9,'p) = Pr {cos 6) e'""^, where Pr{x) is the associated 
Legendre function. This results in certain simplifications, in 
particular, we may only deal with double, not triple, sum- 
mation. We have ^ij^,i2,h=o = 27r$;j,;2, where, according 
to (A5), 



27r 27r 



i,i2 ^ J dwi j dw2x[T{h,h-['\h\),wi\ 



xPi 



cos^o cos^m2 — dSi/dl-i^ 



-i (li W1+I2 W2) 



and ${I,wi,W2) = (27r)-2 J]$,j,,2(-f)e'^'""'+''™''- The 
expression for perturbed DF also retains a mere double sum- 
mation: 

r ^ 1 v^. hdF/dh + hdF/dh i (1^^1+12^2) 

Expression for ^1^,12 i-^) transformed into a more 

compact form. For this purpose we shall first transform 
Pi [cos So cos (102 — dSi/dh)], using the summation theo- 
rem for Legendre polynomials: 

Pi (cos 61 cos 02 — sin 9i sin 62 cos (p) 
I 

= e-"*^P,'=(cos0i)Pr'°(cos6»2). 



Substituting in this formula ^2 = 5 t, 6*1 = | tt — So, 'P ' 
W2 — dS\/dl2 + TT, let us write down 



cos S 



'0 cos^u!2 — dSi/dl2^ 



To close the system, we shall use the integral version of the 
Poisson equation (which appears more convenient for our 
purposes than the Poisson equation itself) 



$(r,S) = 



G 



p{r',e')dV' 



y/j.2 ^ 2 _ 2 rr' cos O 



(^12) 



where dV' is a volume element and O is the angle between 
vectors r and r': cos O — cos S cos S' -l-sin 6 sin 9' cos{(p — (fi'). 
Using (A12) it is simple to obtain the integral connection 
between radial parts of perturbed density p{r) and potential 
x(r). However, the simplest way of obtaining this connection 
is to directly solve, relative to x('')i the known ordinary 
differential equation which follows from the Poisson equation 
after separating the angular dependence : 

1 d r 2 dxjr) ] Ijl + l) , . , n -1 \ I A^^i\ 

Solving Eq. (A13) in terms of Green's function, we shall find 
the required relation: 



47rG 
'2/ + 1 



r' dr p{r') Ti{r,r'). 



(A14) 



where 



ir<y 



r-< =min(r, r'), r-> =max(r, r'). 



For obtaining the integral equation in the desired form it 
is necessary to write down (A14) in action-angle variables 
and to split it into harmonics (hjh)- For this purpose it is 
necessary to select a radial component p{r) from the gen- 
eral expression for density (All) and then to substitute it 
into the r.h.s. of (A14). Further it is necessary to use re- 
lation (AlO), connecting x !i, (2 (-E'l ^) to multiplying 
its both parts by exp [— i {li wi + l2dSi/dl2] and integrat- 
ing over wi. Let us execute the above described procedure. 
We have for p(r): 



J2 Pi'(sinSo)Pr'(0)e-»* 



{ 'W2—dSi/dl2)^ — ik7r 



e-""" (^8) 



Substituting ( A8) into ( A7) and integrating over W2 , we ob- 
tain 



■t>i,j2iE,L) = 2nPl^O)P;'^i8m9o)e"^^Xh,i2iE,L), 



where 



ZTT 

Xi,m{E,L) = j e-'^''^-+'-^^'^"-\[r{E,L,w,)\dw,. 

{AlO) 

For perturbed density p{r, 8, t) — p{r, 6) e we have 
p{r,e)^ f fidv = 



= -iE / dvPl^0)P;'^sm9o)e^''^XHJ2{E,L) 
JidF/dh + I2 dF/dh ^i{iimi+i2VJ2) 

LO — l\Q,\ — l2^2 
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{All) 



p{r') = (Z+ i) / p{r',e')Pi{cose') smO'dO'. (A15) 



Using an explicit form of expression for density (All) and 
substituting the function p{r') found with the help of (A15) 
into r.h.s. of (A14), we obtain: 

x(r)=G^ j dv' j r'^dr' J sinO' dO' Pi {cos 9') J^i{r,r') 



xPf{0)P;'^sm6'o)e^<Xi',i'2{E',L')x 



l[ dF/dlj + 1'2 dF/dl^ 
' Lj-l[ ni{E', L') - l'.^ Q.2{E', L') 



, i{ l'i'w[+l'2'w'2) 



{Aid) 



It is also possible to integrate explicitly over w'2 in (A16). 
For this purpose let us note that in r.h.s. of (A18) there is 
an integration over phase volume J dV' — J dv'dV' , which, 

obviously, may be represented as J dV' — 2n J di'dw'i dw'2. 

Writing down again P;(cosS'), similarly to (A8), as a series. 
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we find 



^7V 

I dw2e"^™^Pi(cose') 



= 2-K P/^ (sin e')P;''^ (0) e""^- e ' ^^i/^-'^ . 

As a result for xC*") we obtain an expression which al- 
ready retains only a single integration over angular variables, 
namely, over w'l : 



2nG J2 I dl'i dl'-, dl's [P'l^ (0) (sin O'^)] 

1 1 , if) 



X 



l[ dF/dlj + 1'2 dF/dl^ 

^ Lj-i[ ni(E', L') - 1'^ n2{E',L') 

X J dw[Ti[r,r'{E',L';w[)] e«{''i 9Si/S^2)_ 



Remembering that F{I) = F{E, L) does not depend on 
Lz ~ I3 we can perform explicitly another integration, 
namely, over 6q. Taking it into account we find 



iv/2 



P''^ (sin e'o)Pi (sin e'o) d (sin Bq) 



-7r/2 



r_i (/ + /2)! 2; + r 

Using a known relation for Pf (0) we obtain for xC**) 



H — — 00 If. — — I 



Xir) 



'1 — ^ '2 



[ fii (E', L')+l'2 Q2 (£', L')] dF/dE' + 1'2 dF/dV 



Lu - /; n^iE', L') - V2 Vl2{E',L') 



X J dw'iJ^i[r,r'{E',L';w'i)] e 


The explicit expression for Z?f is presented in the main text 
(see (|4.4p '). Finally, we obtain the desired set of integral 
equations for Xh,h presented in the main text: 

47rG 

Xh 

zi -1- 1 

Z-, — — 00 i o — — I 



'(p' r/-ii,',\] p»{'i»"i+'2SSi/9-f2) 



4nG ^ >^ I', f dE'LdL' 



xxiii^(i5',i')ni„,,;,;,,^(i5,L;£',L')x 



[z'l ni{E', L') + 1'2 f72 (£',!/')] dF/dE' + 1'2 dF/dL' 
^ Lo-l[ni{E',L')-l'2 Q.2{E\L') ' 

where the kernel is 

n,^.,,^,,,,,(iJ,L;i5',L') = 

= J dwi J dw'^Ti ^r{E,L;'w-i),r'{E',L';w[)\ x 


xexp|i + ;2asi/a/2) - {ii + 1233^/812)]] . 



APPENDIX B. Instabil ity cr iterion for the 
characteristic equation (|4.36p counterpart of 
Penrose - Nyquist plasma criterion 

Let us write down the equation (I4.36|) in the form 
e'-^\z,Q) = 0, where e'''(z,Q) = 1 — Qc{z)/Q and complex 
function Qc{z) is defined as follows: 

dF(x)/dx 



(^) 



dx X ■ 



(Bl) 



(Note that Qc{z) does not depend on Q.) In what follows the 
upper index I is omitted for brevity. It is easy to understand 
that in the points of a real axis z, in which z — Zj = Xj, 
where Xj is any of the extremum points of the function F{x), 
the function Qdz) is real and coincides with the correspond- 
ing value Qc of the "neutral mode". Quotation marks here 
are to reflect the fact that the corresponding value Zj is not 
necessarily a squared frequency of a true neutral mode, since 
the sign of Qc{zj) can be any. Further we shall denote the 
real numbers Qc{zj) as for brevity. 

According to (Bl) the function Qc{z) — Qc{ijJ^), consid- 
ered as a function of ci;, is an analytic function in the tj-plane 
cut along the real axis Im (cu) = 0. From its definition (Bl) 
it also follows that it is continuous and bounded (it tends 
to zero when — ^ oo). We are interested in unstable roots 
of the equations e(cj^, Q) = 0, that is the roots lying in the 
upper half plane uj. The number of such roots, according to 
the argument principle, coincides with the number of poles 
of the function in the upper half plane and is equal 
to A'" = J duje~^de/duj where the contour Cuj is 

shown on 12 (a). Turning to a variable z — ui'^, we find 



N : 



1 

2^ 



de/dz 



dz. 



where the directed contour Cz is the image of a contour 
on a complex plane z. I t is sh own in Fig. 12 (6). 

Following IPenrosd (|l960l ). we shall pass from the com- 
plex plane z to the complex plane Qc- For the number of 
unstable roots N we obtain the following expression 



N = 



27ri 



de dz , 1 
— ■ dz = 

e 2-Ki 



de _ 1 
e 2ni 



dQ, 



(52) 

where the directed contour Cq^ is the image of the directed 
contour Cz on the complex plane Qc- 

Thus the problem is reduced to constructing a contour 
Cq^ . The number of times that this contour encloses in anti- 
clockwise sense the point Qc = Q (lying on the real positive 
half axis of the complex plane Qc) will give us the number 
of unstable roots. We need to formulate rules according to 
which we should image a contour Cz on the plane Qc and to 
establish a direction of motion along it in its various parts. 

1. First of all, note that since the entire circle with a 
big radius on the plane z is imaged into a unique point 
Qc = 0, the entire remaining contour Cq^ on the plane Qc 
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Figure 12. Directed contours of integration Cu and Cz 
plcx planes uj and z respectively. 



corresponds to a horizontal part of the contour Cz , starting 
at a point A, and finishing at a point B. 

2. It is easy to understand that the contour Cq^ is sym- 
metric relative to the horizontal axis Im (Qc) = 0. 

3. Further, we need to understand how the contour Cq^ 
crosses the horizontal axis at a point Qc ~ 0. Note that the 
contour Cq^ should cross the horizontal axis in this point 
(at leasiF^ twice, (i) The first crossing corresponds to the 
start of the contour from Qc = and its finish at Qc = 
when a point z moves from A to B (that is the crossing 
corresponds to coincidence of the "initial" and "final" points 
of the closed contour, (ii) The second crossing corresponds 
to passage of a contour Cz through a point O (that is 2: = 0). 
It follows directly from Equation (Bl), that Qc(0) — 0. Let 
us examine the manner in which the contour Cq^ crosses 
the horizontal axis at Qc = in both cases. The direction - 
up- or downwards - of the crossing of the real axis at this 
point is important here. 

(i) We have from (B2) as \z\ ^ 00 



Qc{z) 



F'{a 



■ dx 



We speak "at least" since it is possible in principle that for 
some Zj corresponding to the extrema of the function F{x), it 
may be true that Qc = 0. 



(x - zy 



■ dx : 



F{x)dx 



M 

z 



M>0. (B3) 



By means of (B3) it is now simple to understand that at a 
point A (that is at z = R — ie, R ^ 00, e — > 0) Im (Qc) — > 
-1-0, and at a point B (that is at 2 = R+ie) Im (Qc) —0. It 
means that this type of crossing is upwards. From (B3) it is 
also visible that crossing occurs so that Re (Qc) +0, that 
is the contour approaches a point Qc = from the right, 
(ti) It is more difficult to establish how crossing occurs that 
corresponds to passage of a point z = on the contour 
Cz. Omitting details, we declare simply that crossing of the 
horizontal axis at a point Qc = in this case also occurs 
upwards - however, the contour approaches this point so that 
Re(Qc) —0, that is from the left. 

4. Now we shall discuss the way the contour crosses 
the horizontal axis at points Qi, corresponding to the ex- 
trema of Fix). Clearly, each such point is crossed twice: the 
first crossing occurs when a point z moves along the pos- 
itive half axis from infinity to the center, and the second 
during subsequent movement in the opposite direction. We 
shall show that both crossings occur in the same direction, 
namely: downwards at points QJ, corresponding to the max- 
ima of F{x), and upwards at points Qc, corresponding to the 
minima. For points z on the real axis z we have 



00 

Qc{z) ^z\^j dx + iTvF'iz) sign (a;)] , 



Im [Qc (z)] — Tvz F' (z) sign (lu) . 



(B4) 



(B5) 

The origin of the factor sign (ui) in (B4) and (B5) may be 
understood from Fig. 12 (b) where it is seen that for points 
z, lying a little below the positive half axis z, i.e., at Im (z) = 
— e, the pole in the integral over x in the right part (B2) is 
indented upwards, and for the points lying a little above the 
real axis, i.e., at lm{z) = e, it is indented downwards. 

Let the point z pass through the corresponding ex- 
tremum point Zj during its first passage, that is from A to 
O along the bottom side of the positive real half axis. Then 
sign (oj) — —1, and the increment Az is negative, Az < 0. 
From (B5) we have 

d 



Aim [Qdz) 



■.F"{z, 



dz 
I Az 



Im [Qdz)] 



Az 
.,)|A2|. 

that is when the 



In its second crossing of this point 
point z moves from O to B along the top side of the positive 
real half axis, we have sign (tj) = -1-1 and Az > 0. So we 
again have 

AIm[Qc(2)] ^ Tvzj F"{zj)Az. 

Thus, indeed, in both crossings the points Qc, corre- 
sponding to extrema, pass in the same manner, namely, the 
minima {F"{zj) > 0) upwards, and the maxima {F"{zj) < 
0) downwards. 

The above recipes are sufficient to construct a directed 
contour Cq^. Having constructed it and found how many 
times it winds around the point Qc ~ Q lying on the hori- 
zontal axis to the right from the origin, we, according to Pen- 
rose's idea, can draw a conclusion on the instability /stability 
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of the system under consideration with a given value of the 
parameter Q. Indeed, from (B2) it is easy to understand that 
if the directed contour Cq^ crosses the horizontal axis from 
below to the right of the point Qc = Q > at least once, 
then instability takes place (A'' ^ 1), since in this case the 
point Qc = Q IS enclosed by the contour in an anticlockwise 
sense. As the parameter Q can have an arbitrary positive 
value, we come to the conclusion that for instability it is 
necessary and sufficient that the contour crosses the hori- 
zontal positive half axis upwards at least once. It means that 
the rightmost crossing of the horizontal axis corresponds to 
a minimum of F{x). This leads us to to formulating the 
instability criterion cited in subsection 4.8 of the main text. 

Note that the task of constructing a contour can be sim- 
plified if we recall that it involves all unstable roots z = zo as 
complex conjugate pairs (this corresponds to pairs of eigen- 
frequencies ujo = ±q + ijS, /3 > 0). Therefore the number 
of times the contour Cq^ winds around a point Qc = Q m 
the complex plane Qc must be even, N — 2L. It is simple 
to understand that L times occur when a point z moves in 
the complex plane z from infinity to zero (along the bottom 
side of the real axis), and L more times occur when it moves 
in the opposite direction (from z = to infinity along the 
top side). Recall that the full contour is symmetric relative 
to the horizontal axis, touches itself at point Qc = 0, and 
the horizontal axis is crossed twice at each point Zj in the 
same direction - either both times downwards (maximum), 
or both times upwards (minimum). Therefore it is enough 
to trace the movement of point z only halfway, say, from 
2: = to infinity and to image only the half of the full con- 
tour which is also a closed contour. For definiteness we shall 
plot only that half which corresponds to movement of point z 
from the center to infinity (along the top side of the real half 
axis). We shall designate this closed directed "semi-contour" 
as Cq^ and, like the full contour, it will start and finish at 
point Qc = 0. Its beginning will be in the second quarter, 
and its end in the fourth quarter of the complex plane Qc- 

Note. The above does not mean that instability exists for any 
< Q < Qc™'"' I where Qc™"'' is a point of right-most crossing 
of a horizontal axis, related to a minimum F(x). If a positive 
Q(^max) g-j^jgj-g g^jgQ £qj. gQjjjg maximum, instability will take place 
only when Q^™''''' < Q < Q^™'""' . Indeed, from Fig. 13 it IS seen 
that in this case the point Qc = Q will only be enclosed if Q lies 
in the specified range. 

This paper has been typeset from a TpjX/ I^T^gX file prepared 
by the author. 
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Figure 13. Schematic sketch of DF F{x) with two maxima 
and one minimum and corresponding arrangement of points Qc 
{j = 1, 2, 3) on axis Im (Qc) = 0. Numbers 0, 1, 2, 3, 4 correspond 
to the order in which the directed semi-contour Cq^ passes over 
corresponding points. The situation when there are two positive 
values Qc with j = 1 and j = 2, corresponding to low maxi- 
mum, and to minimum, is demonstrated. The contour encloses a 
point Qc = Q in anticlockwise sense once provided that the point 
Q is within the range qI^' < Q < Q^^'. These values just are 
boundaries of the range of instability on parameter Q. 
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